We  introduce  a  new  class  of  numerical  differentiation  schemes  constructed  for  the  efficient  solution 
of  time-dependent  PDEs  that  arise  in  wave  phenomena.  The  schemes  are  constructed  via  the 
prolate  spheroidal  wave  functions  (PSWFs).  Compared  to  existing  differentiation  schemes  based  on 
orthogonal  polynomials,  the  new  class  of  differentiation  schemes  requires  fewer  points  per  wavelength 
to  achieve  the  same  accuracy  when  it  is  used  to  approximate  derivatives  of  bandlimited  functions. 
In  addition,  the  resulting  differentiation  matrices  have  spectral  radii  that  grow  asymptotically  as  m 
for  the  case  of  first  derivatives,  and  m2  for  second  derivatives,  with  m  being  the  dimensions  of  the 
matrices.  The  above  results  mean  that  the  new  class  of  differentiation  schemes  is  more  efficient  in  the 
solution  of  time-dependent  PDEs  compared  to  existing  schemes  such  as  the  Chebyshev  collocation 
method.  The  improvements  are  particularly  prominent  in  large-scale  time-dependent  PDEs,  in  which 
the  solutions  contain  large  numbers  of  wavelengths  in  the  computational  domains. 
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1  Introduction 


The  numerical  solution  of  partial  differential  equations  (PDEs)  is  ubiquitous  in  scientific  compu¬ 
tations,  and  as  such  it  is  a  well-developed  subject  that  has  been  widely  studied  (see,  for  example, 
[1,  21,  34,  53,  62]).  One  important  class  of  PDEs  is  the  time-dependent  PDEs,  which  are  PDEs  that 
involve  derivatives  in  both  time  and  spatial  dimensions.  This  kind  of  PDEs  arises  in  the  modeling  of 
physical  phenomena  in  many  scientific  disciplines,  such  as  thermodynamics,  electromagnetics,  and  fluid 
mechanics.  Usually,  the  values  or  the  derivatives  of  the  solution  at  time  t  =  0  (the  initial  conditions) 
and  on  the  boundary  of  the  spatial  domain  (the  boundary  conditions)  are  specified  in  order  to  guarantee 
uniqueness  of  the  solution.  A  common  approach  to  the  solution  of  a  time-dependent  PDE  is  to  first 
approximate  the  spatial  derivatives  by  a  numerical  differentiation  scheme,  discretizing  the  spatial  deriva¬ 
tive  operators.  This  converts  the  PDE  into  a  linear  system  of  ordinary  differential  equations  (ODEs) 
in  time.  The  system  is  then  solved  by  a  numerical  ODE  solver,  such  as  a  Runga-Kutta  scheme  or  a 
predictor-corrector  scheme.  Such  an  approach  for  the  solution  of  time-dependent  PDEs,  which  is  some¬ 
times  referred  to  as  the  “method  of  lines,”  is  studied  in  [54,  55] ,  and  is  applied  in  a  wide  range  of  context 
(see,  for  instance,  [39,  41,  61,  63]). 

The  above  approach  for  the  solution  of  time-dependent  PDEs  requires  the  discretization  of  the  spatial 
derivatives,  for  which  many  numerical  differentiation  methods  are  available.  A  classical  method,  which 
dates  back  to  the  first  half  of  the  nineteenth  century,  is  the  finite  difference  method.  It  approximates 
the  derivative  of  a  function  /  :  R.  — >  R.  at  a  point  x  by  constructing  a  polynomial  p  that  interpolates  /  at 
some  set  of  equispaced  points  around  x,  and  then  taking  the  derivative  of  the  polynomial  p  at  x.  This 
method,  although  simple  to  implement,  has  low  orders  of  accuracy,  and  is  not  suitable  for  the  solution 
of  large-scale  PDEs. 

The  last  few  decades  saw  a  major  advance  in  the  class  of  spectral  methods  for  numerical  differentiation 
and  the  numerical  solution  of  PDEs.  It  was  pioneered  by  the  work  of  Gottlieb,  Orszag,  and  others, 
who  demonstrated  that  spectral  methods  are  particularly  applicable  to  problems  in  fluid  dynamics 
([11,  27,  30]).  Subsequent  development  of  spectral  methods  has  been  contributed  by  the  work  of  Mercier 
[51],  Funaro  [24],  Fornberg  [23],  Trefetlren  [63],  and  Boyd  [8]. 

The  main  idea  of  spectral  methods  is  to  approximate  a  function  f  by  a  finite  series  un  of  smooth 
functions  (f>i , ,(f>n,  and  then  approximate  the  derivatives  of  /  by  the  derivatives  of  the  series  un.  In 
the  context  of  a  time-dependent  PDE,  it  involves  approximating  the  solution  u  of  the  PDE  by  the  finite 
series 

n 

u(x ,  t)  «  un(x,  t)  =  ak(t)(j)k(x),  (1.1) 

k= 1 

where  the  coefficients  a± , . . . ,  an  and  the  functions  ...  ,<f>n  are  chosen  according  to  the  particular  type 
of  spectral  methods  employed.  Among  these  types  of  spectral  methods  are  tau,  Galerkin,  and  collocation 
methods.  The  tau  and  the  Galerkin  methods  involve  converting  the  time-dependent  PDE,  together  with 
its  initial  and  boundary  conditions,  into  a  system  of  ODEs  satisfied  by  the  coefficients  a\ , ...  ,an.  The 
solution  for  aq, . . . , an  then  provides  a  solution  for  u.  When  the  coefficients  a±, . . .  ,an  and  the  functions 
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<f>i , . . . ,  (f>n  are  suitably  chosen,  these  methods  attain  spectral  accuracy,  namely,  the  rates  of  convergence 
of  their  accuracies  depend  only  on  the  smoothness  of  the  solution  u  (see,  for  example,  [30]).  In  particular, 
when  the  solution  u  is  analytic,  the  errors  decay  exponentially  with  the  number  of  functions  n.  Common 
choices  of  (f>i, . . .  ,<j)n  are  trigonometric  functions  for  periodic  problems,  and  Cliebyshev  or  Legendre 
polynomials  for  non-periodic  problems. 

The  collocation  method,  which  is  also  referred  to  as  the  pseudo  spectral  method  (see  [30]),  takes  a 
slightly  different  approach.  First,  a  set  of  points  xi,...,  xm  in  the  spatial  domain,  called  the  collocation 
points,  are  chosen.  Then,  the  coefficients  a\, . . . ,  an  are  computed  such  that  the  series  un  in  (1.1)  exactly 
equals  u  at  the  points  x\, . . .  ,xm,  and  that  un  satisfies  the  boundary  conditions.  This  interpolatory 
procedure  gives  rise  to  an  m  x  m  differentiation  matrix  D  that  takes  the  values  of  u  at  the  collocation 
points  X\, . . . ,  xm  to  the  approximate  values  of  the  spatial  derivatives  of  u  at  aq, . . . ,  xm.  The  “method  of 
lines”  as  described  above  is  then  applied  to  solve  the  time-dependent  PDE,  in  which  the  spatial  derivative 
operators  in  the  PDE  are  discretized  by  these  differentiation  matrices.  Thus,  the  collocation  method 
requires  the  PDE  to  be  satisfied  only  at  the  collocation  points.  For  non-periodic  problems,  the  collocation 
points  are  often  chosen  to  be  the  roots  or  the  extrema  of  the  Chebysliev  or  Legendre  polynomials,  and 
the  interpolants  tfii, . . .  ,<j>n  are  chosen  to  be  polynomials.  For  periodic  problems,  the  collocation  points 
are  often  chosen  to  be  equispaced  points,  and  the  interpolants  are  chosen  to  be  translates  of  the  sine 
function. 

The  collocation  method  amounts  to  the  numerical  approximation  of  spatial  derivatives  of  a  function 
/  at  the  collocation  points  Xi, ,  xrn  via  a  global  interpolant  that  is  exact  at  Xi, ... ,  xm,  and  it  can  be 
viewed  as  high-accuracy  limits  of  finite  difference  methods  (see,  for  instance,  [23]).  Compared  to  the  tau 
and  the  Galerkin  methods,  the  collocation  method  often  leads  to  simpler  systems  of  equations  to  solve, 
while  retaining  spectral  accuracy  (see  [11,  17,  30]);  and  it  is  relatively  easy  to  apply  to  PDEs  that  involve 
variable  coefficients  or  non-linearities.  Therefore,  it  has  become  a  method  of  choice  for  the  numerical 
solution  of  many  types  of  PDEs,  and  by  now  there  is  an  abundance  of  work  in  the  literature  on  both 
its  theoretical  (see,  for  example,  [29,  64,  65,  67,  68])  and  implementational  aspects  (see,  for  example, 
[3,  10,  17,  18,  28,  46]). 

Despite  its  remarkable  accuracy  and  relative  simplicity,  the  collocation  method  has  a  principal  draw¬ 
back  that  limits  its  applicability  in  the  numerical  solution  of  PDEs.  The  differentiation  matrix  D ,  which 
takes  the  values  of  a  function  /  at  the  collocation  points  to  the  approximate  derivatives  of  /  at  the 
same  set  of  points,  is  ill-conditioned  when  non-periodic  boundary  conditions  are  incorporated.  For  the 
case  of  first  derivatives,  the  differentiation  matrix  D  typically  has  a  spectral  radius  of  size  0(N2),  where 
N  is  the  dimension  of  D.  For  the  case  of  second  derivatives,  the  spectral  radius  becomes  0(N4)  (see 
[64,  68]).  This  imposes  strict  stability  requirement  when  the  differentiation  matrix  D  is  combined  with  a 
numerical  ODE  solver  to  solve  a  time-dependent  PDE.  In  many  practical  situations,  the  stability  of  such 
a  combined  scheme  is  determined  by  the  eigenvalues  of  the  differentiation  matrix  D  and  the  time-step 
At  chosen  for  the  ODE  solver.  More  precisely,  the  eigenvalues  of  D ,  multiplied  by  At,  have  to  lie  inside 
the  stability  region  of  the  ODE  solver  in  order  for  the  scheme  to  be  stable.  Therefore,  when  combining 
the  collocation  method  with  an  explicit  time-marching  scheme,  the  time-step  At  typically  has  to  be  of 
size  0(N~2)  when  solving  the  hyperbolic  PDE  =  ux,  and  of  size  0(N~4)  when  solving  the  parabolic 
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PDE  ut  =  uxx.  This  becomes  prohibitive  when  solving  large-scale  PDEs,  for  example,  problems  in 
which  underlying  solutions  span  large  numbers  of  wavelengths  in  the  spatial  domains.  On  the  other 
hand,  using  an  implicit  time-marching  scheme  usually  alleviates  the  time-step  restriction,  but  since  the 
differentiation  matrices  arising  from  the  collocation  method  are  dense,  a  linear  or  non-linear  algebraic 
system  involving  a  dense  matrix  has  to  be  solved  at  each  step,  which  is  expensive.  In  [37],  the  authors 
suggested  a  technique  to  reduce  the  spectral  radii  of  the  first  derivative  matrix  in  the  Chebyshev  colloca¬ 
tion  method  from  0{N2)  to  O(N)  by  first  transforming  the  collocation  points.  The  technique  has  been 
applied  in  a  number  of  situations  (see,  for  example,  [4,  18,  32,  50]);  however,  it  requires  careful  choice 
of  a  transformation  parameter  in  order  to  maintain  desired  accuracy.  Moreover,  as  pointed  out  in  [50], 
for  practical  values  of  N,  the  anticipated  0(N~1)  time-step  size  is  not  attained,  and  is  replaced  by  at 
most  doubling  of  the  time-step  allowed  by  the  Chebyshev  collocation  method.  Also,  the  transformation 
destroys  the  quadrature  properties  of  the  roots  of  the  Chebyshev  polynomials,  which  is  undesirable  in 
certain  applications  (see  [13]). 

Recently,  there  has  been  a  growing  interest  in  the  development  of  spectral  and  pseudospectral  methods 
based  on  the  prolate  spheroidal  wave  functions  (PSWFs),  which  were  introduced  in  a  series  of  work  by 
Slepian  et.  al  in  the  context  of  the  analysis  of  bandlimited  functions  defined  on  intervals  (see  [42,  43,  57, 
58,  59]).  In  [69],  the  authors  constructed  quadrature  and  interpolation  formulas  for  bandlimited  functions 
based  on  the  PSWFs,  demonstrating  that  the  PSWFs  are  a  natural  tool  for  the  design  of  numerical 
algorithms  for  bandlimited  functions.  In  particular,  the  PSWFs  and  their  variants  have  been  used  as 
basis  functions  in  the  construction  of  spectral  and  pseudospectral  methods  for  the  solution  of  ODEs  and 
PDEs  describing  wave  phenomena  (see,  for  example,  [6,  9,  13,  38,  39,  45,  66]).  For  instance,  the  authors 
in  [38,  45]  constructed  collocation  methods  with  the  PSWFs  as  basis  functions  to  solve  the  Schrodinger’s 
equation,  and  the  authors  in  [13]  constructed  collocation  methods  based  on  quadrature  nodes  and  roots 
associated  with  the  PSWFs,  and  apply  them  to  the  solution  of  hyperbolic  PDEs.  More  recently,  the 
authors  in  [6]  developed  a  two-dimensional  solver  for  wave  equations  that  involves  spatial  derivative 
operators  constructed  using  bases  of  “approximate”  PSWFs.  All  of  the  above  results  reinforced  the 
observations  in  [69]  that  for  problems  involving  bandlimited  functions,  collocation  methods  based  on  the 
PSWFs  require  fewer  points  per  wavelength  to  achieve  the  same  accuracy  when  compared  to  collocation 
methods  based  on  orthogonal  polynomials,  such  as  Chebyshev  or  Legendre  polynomials.  However,  while 
these  results  are  encouraging,  the  development  of  pseudospectral  methods  based  on  the  PSWFs  is  still  in 
its  nascent  stage.  First,  further  investigations  into  the  relationship  between  the  choice  of  the  PSWFs  in 
the  approximation  (1.1)  and  the  resulting  accuracies  in  numerical  differentiation  and  solution  of  PDEs 
are  needed.  Second,  while  preliminary  results  in  [9,  13,  39,  69]  indicate  that  collocation  methods  based 
on  PSWFs  lead  to  differentiation  matrices  of  somewhat  smaller  condition  numbers  compared  to  those 
based  on  Chebyshev  or  Legendre  polynomials,  a  systematic  picture  on  the  condition  numbers  of  the 
differentiation  matrices  in  relation  to  the  choice  of  the  PSWFs  and  the  accuracy  requirement  is  still 
lacking. 

In  this  paper,  we  introduce  a  new  class  of  numerical  differentiation  schemes  constructed  using  the 
PSWFs.  The  schemes  are  constructed  based  on  the  approximation  of  a  function  /  by  a  linear  combination 
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of  PSWFs: 


n—  1 

f(x)  ~'52aji’j(x),  (1-2) 

i= o 

and  differ  from  existing  schemes  based  on  PSWFs  in  two  main  aspects.  First,  the  PSWFs  in  the  approx¬ 
imation  (1.2)  are  determined  by  the  choice  of  the  bandlimit  parameter  c  and  an  accuracy  requirement  e. 
In  particular,  the  number  of  PSWFs  n  in  (1.2)  is  determined  by  c  and  e.  This  is  different  from  schemes 
in  [9,  13,  38,  45],  in  which  the  bandlimit  parameter  c  and  the  number  of  functions  n  are  independently 
(and,  at  times,  somewhat  arbitrarily)  chosen.  The  second  difference  is  that,  as  opposed  to  existing 
pseudospectral  methods,  which  often  force  the  interpolant  of  the  function  /  to  exactly  equal  to  /  at  a 
set  of  collocation  points,  we  do  not  force  the  approximation  (1.2)  to  exactly  hold  at  a  set  of  points  in 
the  computation  of  the  coefficients  a.\ ,...,an.  Instead,  a  least-squares  type  procedure,  based  on  the 
quadratures  for  bandlimited  functions  constructed  in  [69],  is  used  to  compute  the  coefficients  a\ , ,an. 
The  end  result  is  an  m  x  m  differentiation  matrix  that  takes  the  values  of  the  function  /  on  a  chosen  set 
of  quadrature  nodes  xm  to  the  approximate  values  of  the  derivatives  of  /  on  x±, . . . ,  xm.  The  only 

requirement  on  the  quadrature  is  that  it  integrates  the  products  of  the  PSWFs  ipQ, . . . ,  ^-1  to  sufficient 
accuracy.  In  particular,  we  do  not  require  the  number  of  PSWFs  n  in  (1.2)  to  equal  the  number  of  nodes 
m  used  in  the  construction  of  the  differentiation  matrix,  or  prescribe  any  other  functional  relationship 
between  n  and  in. 

One  benefit  of  the  new  class  of  numerical  differentiation  schemes,  which  was  indicated  in  the  discussion 
in  [69],  is  that  when  dealing  with  problems  involving  bandlimited  functions,  it  requires  fewer  points  per 
wavelength  to  achieve  a  prescribed  accuracy,  compared  to  schemes  based  on  orthogonal  polynomials, 
such  as  the  Chebysliev  collocation  method.  More  importantly,  in  the  solution  of  time-dependent  PDEs 
with  non-periodic  boundary  conditions,  the  resulting  first  and  second  derivative  matrices  have  spectral 
radii  that  grow  as  c  and  c2  respectively,  for  fixed  e.  This  corresponds  to  differentiation  matrices  with 
spectral  radii  that  grow  as  m  and  m2  respectively  for  large  m,  with  m  being  the  dimensions  of  the 
matrices.  The  result  means  that,  when  we  combine  these  differentiation  matrices  with  an  explicit  time¬ 
marching  scheme  to  solve  the  PDE,  a  larger  time-step  At  can  be  chosen  compared  to  when  the  collocation 
method  is  used,  while  maintaining  stability.  Combining  with  the  result  that  differentiation  matrices  of 
smaller  dimensions  are  needed  to  achieve  the  same  accuracy  compared  to  the  collocation  method,  the 
new  class  of  differentiation  schemes  are  more  efficient  in  the  solution  of  time-dependent  PDEs  involving 
bandlimited  functions.  The  advantage  is  critical  in  the  case  of  large-scale  PDEs,  in  which  a  large  number 
of  points  are  needed  to  discretize  the  solution  in  the  spatial  domain. 

The  paper  is  structured  as  follows.  Section  2  provides  the  necessary  mathematical  and  numerical 
preliminaries.  Section  3  describes  the  construction  of  the  new  class  of  numerical  differentiation  schemes 
based  on  the  PSWFs,  as  well  as  modifications  to  the  schemes  when  boundary  conditions  are  incorporated. 
In  Section  4,  we  present  numerical  results  pertaining  to  the  accuracy  and  stability  properties  of  the 
schemes,  and  discuss  the  results  of  several  numerical  experiments  when  the  schemes  are  applied  to  the 
solution  of  time-dependent  PDEs  and  the  associated  eigenvalue  problems.  Finally,  Section  5  summarizes 
the  work  and  discusses  possible  extensions  and  generalizations. 
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2  Mathematical  and  numerical  preliminaries 

2.1  Quadrature  and  interpolation 

2.1.1  Generalized  Gaussian  quadratures 

The  quadrature  rules  considered  in  this  paper  are  of  the  form 

n 

(2-1) 

3=1 

where  the  points  Xj  £  R.  and  coefficients  Wj  £  M  are  referred  to  as  the  nodes  and  weights  of  the 
quadrature,  respectively.  They  serve  as  approximation  to  integrals  of  the  form 

f  (f>(x)oj(x)dx,  (2-2) 

J  a 

where  w  :  [a,  b]  — >  R.  is  an  integrable  non-negative  function. 

Quadratures  are  typically  chosen  so  that  the  quadrature  (2.1)  is  equal  to  the  integral  (2.2)  for 
some  set  of  functions,  commonly  polynomials  of  some  fixed  order.  One  main  example  is  the  classical 
Gaussian  quadrature,  which  consists  of  n  nodes  and  integrates  polynomials  of  degree  up  to  2 n  —  1 
exactly.  The  notion  of  Gaussian  quadratures  can  be  generalized  to  other  systems  of  functions  as  follows 
(see  [14,  47,  70]): 

Definition  2.1.  A  Gaussian  quadrature  for  a  set  of  2 n  functions  <j>i, . . . ,  </>2„  :  [a,  b]  — >■  R.  with  respect  to 
a  weight  function  co  :  [a,  6]  — >  R+  is  a  quadrature  rule  with  n  weights  and  nodes  that  integrates  exactly 
(f>i  with  respect  to  w  for  all  i  =  1, . . . ,  2 n.  The  weights  and  nodes  of  a  Gaussian  quadrature  will  be 
referred  to  as  Gaussian  weights  and  nodes,  respectively. 

Remark  2.1.  While  the  existence  of  generalized  Gaussian  quadratures  has  been  proven  for  a  fairly  broad 
class  of  systems  of  functions  for  more  than  100  years  (see,  for  instance,  [36,  40,  48,  49]),  the  constructions 
found  in  [25,  35,  36,  40,  48,  49]  do  not  easily  yield  numerical  algorithms  for  the  design  of  such  quadratures. 
Such  algorithms  have  been  constructed  recently  (see  [14,  47,  70]). 


2.1.2  Discretization  of  square  integrable  functions 

We  shall  say  that  a  quadrature  rule  with  nodes  xi,...,xn  £  [a,  6]  and  positive  weights  Wi, . . .  ,wn 
discretizes  a  collection  of  square  integrable  functions  /i, . .  • ,  fm  defined  on  the  interval  [a,  b]  if  it  integrates 
exactly  all  pairwise  products  of  /i, . . . ,  /m;  in  other  words,  if 

/i b  n 

fi(x)fj(x)dx  =  fi{xi)fj(xi)wi  (2.3) 

;=i 


holds  for  all*,  j  =  1, . . . ,  to. 
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If  Xi, . . . ,  xn,  w i, . . . ,  vun  is  a  quadrature  discretizing  a  collection  of  functions  /i, . . . ,  fm  in  L2([a,  &]), 
then  the  map  T  from  the  span  S  of  /i , . . . ,  fm  to  the  Euclidean  space  ln  taking  the  function  f  to  the 
vector 


/  /(*i)v^T  \ 

f(x2)y/w^ 

V  f{xn)y/Wn  / 


(2.4) 


is  a  Hilbert  space  isomorphism  (a  bijection  which  preserves  inner  products)  of  the  subspace  S  onto  the 
subspace  of  the  Euclidean  space  R"  spanned  by  the  vectors 


I  \ 

/  fm{x\)yfw{  \ 

fl(%2)y/W2 

,  .  .  .  , 

fm{X2)y/u^ 

V  h{xn)y/Wn  ) 

V  fm (xn)y/ Wrl  J 

(2.5) 


2.1.3  Stable  interpolation  on  quadrature  nodes 

If  u\, . . . ,  Uk  is  a  collection  of  orthonormal  functions  in  E2([a,  6]),  and  x\, . . . ,  xn,  Wi, . . .  wn  is  a  quadrature 
discretizing  w1; . . . , Uk,  then  Xi, ...  ,xn  serve  as  stable  interpolation  nodes  for  the  span  of  u±, ...  ,Uk-  In 
particular,  for  a  function  /  defined  on  [a,  6],  we  can  compute  stably  the  coefficients  a±  ,...,otk  in  the 
linear  combination 

k 

f  =  Y]  otjUj  (2.6) 

using  the  values  /i ,  . . . ,  fn  of  /  at  the  nodes  x\, ...  ,xn.  Let  U  be  the  n  x  k  matrix  with  entries 


Ui.j  —  Uj  (^i )  i 

and  let 

F=  (fl,  ■  ■■,  fn)T , 

a  =  (ai, . . .  ,ak)T, 

then  (2.6)  implies  that 

F  =  Ua. 

Multiplying  both  sides  of  (2.10)  by  the  n  x  n  diagonal  matrix  W  with  entries 

Witi  =  y/m,  i  =  1, . . .  ,n, 


(2.7) 

(2.8) 
(2.9) 

(2.10) 

(2.11) 


we  obtain 

WF  =  Aa,  (2.12) 
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where 


A  =  WU.  (2.13) 

Since  the  quadrature  aq, . . . ,  xn,  uq, . . . ,  wn  discretizes  u±, . . . ,  «*,,  the  matrix  A  has  orthonormal  columns. 
Therefore,  (2.12)  provides  a  numerically  stable  formula  to  compute  the  coefficients  a: 


a  =  A*WF.  (2.14) 

In  other  words,  a  can  be  obtained  by  applying  to  F  a  diagonal  matrix  followed  by  a  matrix  with 
orthonormal  rows.  Values  of  the  function  /  at  any  point  x  €  [a,b\  can  then  be  computed  using  (2.6), 
provided  that  a  scheme  for  evaluating  the  functions  U\, . . .  ,Uk  at  arbitrary  points  is  available. 

Remark  2.2.  For  simplicity  of  discussion,  we  have  assumed  that  the  representation  (2.6)  of  /  holds 
exactly.  In  actual  computations,  we  often  represent  /  by  a  truncated  series  of  the  form  (2.6)  that  is 
accurate  to  a  precision  e.  As  long  as  the  quadrature  rule  x\, . . . ,  xn,  w\, . . . ,  wn  discretizes  U\ , . . . ,  Uk, 
(2.14)  is  a  stable  interpolation  formula  for  /  that  is  accurate  to  e. 

Remark  2.3.  At  first  glance,  the  above  procedure  for  the  computation  of  a  seems  rather  limited  in  scope, 
since  it  relies  on  the  quadrature  being  exact  for  all  pairwise  products  of  u\, . . . ,  Uk-  However,  as  long  as 
the  quadrature  is  reasonably  accurate  for  all  pairwise  products  of  U\, . . .  ,Uk ,  the  matrix  A  in  (2.12)  is 
close  to  having  orthonormal  columns,  and  therefore  sufficiently  well-conditioned.  In  this  case,  a  can  be 
computed  stably  by  solving  (2.12)  using  the  least-squares  method. 


2.2  Prolate  spheroidal  wave  functions 

2.2.1  Basic  facts 

In  this  subsection,  we  summarize  some  basic  facts  about  PSWFs.  Unless  stated  otherwise,  all  of  these 
facts  can  be  found  in  [42,  59,  69]. 

Given  a  real  c  >  0,  we  denote  by  Fc  the  operator  L2([— 1, 1])  — >  L2([—  1, 1])  defined  by  the  formula 

Fc{(t>){x )  =  J  eicxt<j)(t)dt.  (2.15) 

Obviously,  Fc  is  compact.  We  denote  the  eigenvalues  of  Fc  by  A0,  Ai, . . . ,  A  j, . . .  such  that  |  Ay_i  |  >  |  Ay  | 

for  all  integer  j  >  1.  For  each  integer  j  >  0,  we  denote  by  ipj  the  eigenfunction  corresponding  to  Ay ;  in 

other  words,  the  integral  equation 

Xjtpj(x)  =  J  elcxtipj(t)dt  (2.16) 

holds  for  all  x  G  [—1,1].  Following  [69],  we  adopt  the  convention  that  -i/’j  are  normalized  such  that 

||VblU2([-MD  =  1  for  a11  J- 

The  following  theorem  summarizes  the  basic  properties  of  t/y  and  A  j ,  and  can  be  found  in  a  slightly 
different  form  in  [69]. 


Theorem  2.1.  For  any  real  c  >  0,  the  eigenfunctions  ifo,  ifh  ■  ■  ■  >  °f  the  operator  Fc  are  purely  real, 
orthonormal ,  and  complete  in  i2  ( [ —  1 , 1] ) .  V’j  is  even  for  all  even  j  and  is  odd  for  all  odd  j.  Each 
function  ipj  has  exactly  j  simple  roots  in  (—1,1).  All  eigenvalues  Xj  of  Fc  are  non-zero  and  simple;  X j 
is  purely  real  for  all  even  j  and  is  purely  imaginary  for  all  odd  j;  in  particular,  Xj  =  *J'|Aj|. 

We  define  the  self-adjoint  operator  Qc  :  L2{{—  1, 1])  — >  L2{[—  1, 1])  by  the  formula 

QM)  =  -  I'  sm(c(a;-f))0(t)df.  (2.17) 

7Tj_i  X~t 


A  simple  calculation  shows  that 

Qc=^'K-  Fc,  (2.18) 

and  that  Qc  has  the  same  eigenfunctions  as  Fc.  Moreover,  the  jtli  (in  descending  order)  eigenvalue  /ij 
of  Qc  is  related  to  A  j  by  the  formula 

W  =  ^-|Ai|2.  (2.19) 

The  operator  Qc  is  closely  related  to  the  operator  Pc  :  L2(R)  — >  L2(K)  defined  by  the  formula 


1 

7T 


sin(c(x  —  t)) 
x  —  t 


4>(t)dt, 


(2.20) 


which  is  the  well-known  orthogonal  projection  operator  onto  the  space  of  functions  of  band  limit  c  defined 
on  R. 

The  following  theorem  describes  the  behavior  of  the  spectrum  of  Qc,  and  can  be  found  in  [69].  It  is 
proven  in  a  slightly  different  form  in  [44]. 

Theorem  2.2.  For  any  c  >  0  and  0  <  a  <  1  the  number  N  of  the  eigenvalues  Pj  that  are  greater  than 
a  satisfies  the  equation 

N  =  —  +  (- 2  log  - — log(c)  +  0(log(c)).  (2.21) 

7T  y  71  a  J 

Equation  (2.21)  implies  that  for  large  c  >  0,  Qc  has  about  2c/-7T  eigenvalues  with  magnitudes  very 
close  to  1,  followed  by  order  log(c)  eigenvalues  decaying  exponentially  from  1  to  nearly  0;  the  rest  of  the 
eigenvalues  are  all  very  close  to  0. 

The  eigenfunctions  ipo,ipi, . . .  ,ipj, . . .  of  Qc  turn  out  to  be  the  PSWFs,  a  fact  well  known  from 
classical  mathematical  physics  (see,  for  example,  [52]).  The  following  theorem  formalizes  this  statement, 
and  is  proven  in  a  more  general  form  in  [57]. 

Theorem  2.3.  For  any  c  >  0,  there  exists  a  strictly  increasing  unbounded  sequence  of  positive  numbers 
Xo>  Xii  •  •  •  such  that,  for  each  integer  j  >  0,  the  differential  equation 


(1  —  x2)ip  ( x )  —  2xif  (. x )  +  (xj  —  c2x2)if>(x)  =  0 


(2.22) 


has  a  solution  that  is  continuous  on  [—1,1].  Moreover,  all  such  solutions  are  constant  multiples  of  ipj 
(2.16). 
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Remark  2.4.  To  be  more  precise,  the  eigenfunctions  ipo,ipi,  ■  ■  ■  corresponding  to  the  operator  Fc  (and 
Qc)  should  be  denoted  ipQ,ip(, ....  For  simplicity  of  notation,  we  will  omit  the  superscript  c  wherever 
the  value  of  c  is  clear  from  the  context. 


2.2.2  Approximations  of  bandlimited  functions  by  PSWFs 

In  this  subsection,  we  define  bandlimited  functions  and  provide  a  brief  review  on  their  approximation 
by  PSWFs.  More  details  can  be  found  in  [56,  69]. 

Definition  2.2.  Let  I  be  an  interval  in  the  real  line.  A  function  /  :  I  — »•  R  is  said  to  be  bandlimited  if 
there  exists  a  positive  real  c  and  a  function  a  £  L2([—  1, 1])  such  that 

/(*)  =  J  eicxta(t)dt  (2.23) 

for  all  x  £  I.  Moreover,  in  this  case  /  is  said  to  have  bandlimit  c. 

Obviously,  (2.16)  implies  that  ipj  is  bandlimited  by  c  for  all  integer  j  >  0. 

Since  the  PSWFs  ip o,ip(,  ■  ■  ■  constitute  a  complete  orthonormal  basis  in  L2{{—  1, 1])  (see  Theorem 
2.1),  a  function  /  :  [— 1, 1]  — »•  R  with  bandlimit  c  has  the  infinite  series  expansion 


OO 


f(x)  =]Ta^0r), 

3=0 


(2.24) 


where  the  coefficients  ay  are  given  by  the  formula 


oij  =  J  ipj(x)f(x)dx, 


(2.25) 


for  all  integer  j  >  0.  The  series  (2.24)  is  convergent  in  L  ,  moreover,  it  is  shown  in  [56]  that  the  series  in 
(2.24)  converges  uniformly  to  /  on  [—1, 1],  The  following  lemma,  which  is  proven  in  a  slightly  different 
form  in  [56],  provides  a  bound  for  approximating  a  bandlimited  function  with  a  truncated  expansion  in 
PSWFs. 


Lemma  2.4.  Suppose  that  f  :  [—1,1]  — >  R  has  bandlimit  c  and  is  expressible  as  (2.23)  and  (2.2)). 
Then,  for  any  non-negative  integer  k,  the  bound 


k 

3=0 


(2.26) 


holds  for  all  x  £  [—1, 1],  where 


£c=  max  \tpj(x)\,  j  =  1,2,.... 

J  — 1<X<1  J 


(2.27) 


The  bound  (2.26)  indicates  that  the  error  is  roughly  of  the  order  of  |Afc|,  provided  that  k  is  in  the 
range  where  the  eigenvalues  of  Qc  are  decaying  exponentially  (see  Theorem  2.2).  The  numerical  results 
in  [69]  confirm  this  observation. 

2.2.3  Quadratures  for  bandlimited  functions 

We  shall  say  that  a  quadrature  rule  with  nodes  x±, . . .  ,xn  €  [—1,1]  and  positive  weights  wi, . . .  ,wn 
integrate  functions  on  [— 1, 1]  with  bandlimit  c  to  precision  e  if  for  any  function  /  on  [—1,1]  of  the  form 
(2.23),  we  have 

1  n 

f{x)dx-J2f(xj)wj 

_1  i= i 

One  procedure  to  construct  Gaussian  quadratures  for  bandlimited  functions  is  to  use  the  Newton-type 
nonlinear  optimization  algorithm  of  [14].  Specifically,  for  bandlimit  c  and  precision  e,  the  algorithm 
constructs  an  n/ 2-point  Gaussian  quadrature  integrating  exactly  the  first  n  PSWFs  i />§,  V’l,  •  •  ■ ,  V’n-i, 
where  n  is  the  smallest  integer  such  that  the  corresponding  eigenvalue  has  magnitude  less  than  e,  i.e. , 

|A„_i|  >e>  |A„|.  (2.29) 

The  constructed  quadrature  then  integrates  all  functions  on  [—1,1]  with  bandlimit  c  to  precision  e, 
provided  n  is  in  the  range  where  the  eigenvalues  of  Qc  decay  exponentially. 

Remark  2.5.  The  procedure  described  above  is  expensive,  and  requires  order  n3  operations  (see  [14]). 
An  alternative  procedure  is  described  in  [69],  and  is  based  on  a  generalization  of  the  Euclid’s  division 
algorithm  to  PSWFs.  While  it  is  less  expensive  than  the  procedure  of  [14],  the  quadrature  constructed 
is  a  bit  less  accurate.  Numerical  examples  demonstrating  the  performance  of  quadratures  constructed 
by  both  procedures  can  be  found  in  [69]. 

2.2.4  Interpolation  of  bandlimited  functions 

As  discussed  in  Section  2.2.2,  given  a  function  /  :  [—1, 1]  — >  R.  with  bandlimit  c,  we  can  approximate  it 
by  a  linear  combination  of  the  first  n  PSWFs  ip^,  ■  •  • ,  ip n-i '■ 

n—  1 

Eai^.  (2-3°) 

3=  0 

with  an  error  roughly  of  the  order  of  |A„|,  provided  that  n  is  in  the  range  where  the  eigenvalues  of 
Qc  decay  exponentially.  Given  a  quadrature  xi,...  . . .  ,wm  that  integrates  exactly  all  pairwise 

products  of  ip g, . . . ,  V’n-i;  we  can  apply  same  argument  as  in  Section  2.1.3  to  express  the  problem  of 
computing  the  coefficients  ay  as 

WF  =  Aa,  (2.31) 


<£J  \cr(t)\dt. 


(2.28) 
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where 


(2.32) 

(2.33) 


F  =  {f{x  1),  •  •  ■J(xm))T, 

OL  =  («o,  •  •  •  ,  0Ln- l)T, 

W  is  an  m  x  to  diagonal  matrix  with  entries 

Witi  =  i  =  1, . . . ,  to,  (2.34) 

and  ^4  is  an  to  x  n  matrix  with  entries 

A,j  =  (2.35) 

The  matrix  A  has  orthonormal  columns,  and  a  can  be  computed  stably  using  the  formula 


a  =  A*WF. 


(2.36) 


(2.30)  then  gives  an  interpolation  formula  for  /  with  error  proportional  to  |A„|. 

Remark  2.6.  Theorem  7.1  of  [69]  shows  that  a  quadrature  integrating  all  functions  on  [—1, 1]  of  bandlimit 
2c  to  precision  |A°|2  guarantees  that  the  matrix  A  in  equation  (2.31)  is  close  to  having  orthonormal 
columns,  and  hence  sufficiently  well-conditioned  for  the  accurate  computation  of  a  (see  Remark  2.3). 
In  the  numerical  experiments  in  Section  4.1.1,  we  use  a  quadrature  corresponding  to  bandlimit  2c  and 
precision  |A“  |2  x  10-1°  in  order  to  ensure  that  A  has  orthonormal  columns. 

Remark  2.7.  Both  the  procedures  of  [14]  and  [69]  (see  Remark  2.5)  can  be  used  to  construct  quadratures 
for  the  interpolation  of  bandlimited  functions.  In  addition,  the  numerical  results  in  [69]  show  that 
the  two  quadratures  have  virtually  identical  performance  when  used  for  interpolation  (as  opposed  to 
integration),  achieving  the  same  accuracy  with  the  same  number  of  nodes. 

2.3  Numerical  differentiation 

2.3.1  Finite  difference  method 

The  finite  difference  method  is  a  classical  numerical  differentiation  method  based  on  local  polynomial 
interpolation  of  a  function  at  equispaced  points.  Given  a  function  /  :  K  — >  R  and  li  /  0,  the  Taylor 
expansion  of  f(x  +  h)  around  the  point  x: 

f(x  +  h)  =  f(x)  +  +  0(h2) 

gives  rise  to  the  following  one-sided  difference  formula  for  the  approximation  of  fix) 

rr  f{x  +  h)-f{x ) 


(2.37) 


(2.38) 
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which  has  an  error  of  0{h).  Note  that  formula  (2.38)  also  results  from  linear  interpolation  of  /  at  the 
points  x  and  x  +  h.  By  constructing  the  degree-two  polynomial  p  that  interpolates  /  at  the  points  x  —  h, 
x,  and  x  +  h,  and  then  computing  p'(x),  we  obtain  the  centered  difference  formula  for  f'(x): 


rr  fix  +  h)  -  f(x  -  h) 

=  - TyT - 


(2.39) 


which  has  an  error  of  0(h2);  in  other  words,  (2.39)  is  second-order  accurate.  Alternatively,  (2.39)  is 
derived  by  combining  the  Taylor  expansions  of  f(x  +  h)  and  f(x  —  h)  around  the  point  x. 

Applying  the  same  idea  of  local  polynomial  interpolation,  we  obtain  finite  difference  formulas  with 
higher  orders  of  accuracy,  or  for  higher-order  derivatives.  Table  1  lists,  with  corresponding  orders  of 
accuracy,  some  finite  difference  formulas  for  the  approximations  of  f  and  f"  using  the  values  of  /  at 
equispaced  points  of  step-size  h.  More  comprehensive  tables  of  finite  difference  formulas  can  be  found, 
for  instance,  in  [22,  23]. 


Table  1:  Finite  difference  formulas  for  the  approximations  of  f  and  f"  using  the  values  of  /  at  equispaced 
points. 


Derivative 

/'(*) 

/'(*) 

/'(*) 

/'(*) 

/'(*) 

/"(*) 

/"(*) 

/"(*) 

/"(*) 

/"(*) 


Finite  difference  formula 


f(x+h)-f(x-h) 

2  h 

-3f(x)+4f(x+h)-f(x+2h) 

2  h 

f(x-2h)  —  8f{x—h)-\-8f(x-\-h)—f(x-\-2h) 

12 h 

—3/ (x—h)  —  10f  (x) 

12  h 

— 25/(x)+48/(x+/i)— 36  f  (x+2h)-\-16  f  (x-\-3h) — 3f(x-\-4h) 
12/?, 


v)+18  fix +h) — 6  f  (x-\-2h)-\-  f  (x+3h) 


f(x-h)-2f(x)+f(x+h) 
h 2 

2f(x)—5f(x+h)+4:f(x+2h)  —  f(x+3h) 
h 2 

-f(x—2h)+16f(x—h)—30f(x)+16f(x-\-h)—f(x+2h) 

12h? 

10f(x—h)  —  15f(x)—4f(x+h)+14f(x-\-2h)—6f(x+3h)+f(x+4h) 

12 h? 

45f(x)-154f(x-\-h)-\-214f(x-\-2h)  —  156f(x-\-3h)-\-61f(x-\-4h)  —  10f(x-\-5h) 

12 h? 


Order  of  accuracy 
2 
2 
4 
4 
4 
2 
2 
4 
4 
4 


2.3.2  Chebyshev  collocation  method 

In  this  subsection,  we  briefly  outline  the  Chebyshev  collocation  method  for  numerical  differentiation. 
Detailed  discussions  of  its  implementation,  as  well  as  its  accuracy  and  stability  properties,  can  be  found 
in,  for  instance,  [8,  23,  30,  63]. 

Consider  the  Chebyshev-Gauss-Lobatto  collocation  points  x\, . . . , xrn  defined  on  [—1,1]  via  the  for¬ 
mula 

Xi  =  cos 

In  particular,  xi,...,xm  are  arranged  in  descending  order,  and  we  have  x,\  =  1  and  xm  =  —  1.  The 


1) 

m  —  1 


,  i  =  1, . . . ,  to. 


(2.40) 
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points  X\, . . . ,  xm  are  the  extrema  of  the  Chebyshev  polynomial  of  order  to  —  1: 

Tm_i(x)  =  cos  ((m  —  1)  cos-1  x)  . 


(2.41) 


Given  a  function  /  :  [—1,1]  — >  R,  the  Chebyshev  collocation  method  approximates  the  derivatives  of 
/  at  the  collocation  points  xi,...  ,xm  by  first  constructing  the  interpolating  polynomial  fm  of  degree 
to  —  1  that  agrees  with  f  at  Xi, ... ,  xm,  i.e. , 

m 

fm(x)  =  f{x3)gj{x),  (2.42) 

i=i 

where  g3  is  a  polynomial  of  degree  to  —  1  uniquely  defined  by  its  values  at  xi, . . . ,  xm: 


9j{xk) 


1  if  j  =  k, 
0  if  j  ^  k. 


(2.43) 


It  can  be  shown  (see,  for  example,  [8])  that 


,  ,  (-1  Vil-x^T^x)  .  i 

9j{ x)  = - 7 - tuu - ,  J  =  l,...,m, 

Cj{m-  1  y{x-Xj) 


where 


c,  = 


2  if  j  =  1,  to, 

1  if  j  =  2, . . . ,  to  —  1. 


(2.44) 


(2.45) 


The  interpolating  polynomial  (2.42)  provides  approximations  d\, . . . ,  drn  of  the  first  derivatives  of  /  at 
Xi, ... ,  xm  by  the  formula 

m 

di  =  '^2f(xj)g'j(xi),  i  —  1, . . . ,  m.  (2.46) 

j= i 

Therefore,  the  to  x  to  Chebyshev  differentiation  matrix  D  with  entries  D,  j  =  g'j(xi)  is  a  linear  operator 
taking  the  values  of  the  function  /  at  the  collocation  points  x±, . . . ,  xm  to  the  approximate 

values  di, ... ,  dm  of  f  at  X\, . . . .  xm.  A  straightforward  calculation  gives  us  an  explicit  formula  for  the 
entries  of  D: 


Du 

Di,i 


1  Xi 


2  1  -  xf 
2  (to  —  l)2  +  1 
6 


,  i  ±  1,  to, 

=  -Dr, 


(2.47) 

(2.48) 

(2.49) 


Alternatively,  since  the  interpolating  polynomial  (2.42)  to  the  vector  (1, . . . ,  1)T  is  exactly  the  constant 
function  /(x)  =  1,  and  that  f'(x )  =  0  for  all  x,  the  matrix  D  must  map  (1, . . . ,  1)T  to  the  zero  vector, 
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by  construction.  Thus,  we  can  compute  the  diagonal  entries  of  D  by  the  formula 


m 

33 1 .  /  ^  ^  i  1 5  -  *  -  5  ^  C 


(2.50) 


The  reduction  of  rounding  errors  in  computing  the  diagonal  entries  of  D  via  (2.50)  instead  of  (2.48)-(2.49) 
are  discussed  in  [2,  3]. 

By  the  same  argument,  the  m  x  to  Chebysliev  differentiation  matrix  D  for  the  approximation  of 
second  derivatives  of  functions  /  :  [—  1, 1]  — >  R.  has  entries  g"(xi).  The  off-diagonal  entries  of  D  are  given 
by  the  formula 


fS  1  /  ,  \  j-i-,'-)-!  2  XiXi  Xj  ...  -  /  i 

Di,j  —  ~(  —  1)  7i  _.2\  / ITT?’  1  7^  J'l  1  7^  1)TO) 


(1  ~  xi)(xi  —  xj)2 


Dl,j  - 


Dm,j  — 


(_1)J  /— 4(m  —  l)2  -  2 


+ 


Cj  V  3(^1  -  x3 )  1  Oi  -  Zj)2  '  ’  J  ^ 

(-l)J  /(-l)m"1(4(TO- l)2 +  2)  (-l)™-^ 


3(a;m  -  Zj)  (xm  -  a^)2 

The  diagonal  entries  of  D  are  then  computed  by  the  formula 


,  3  ±  rn. 


(2.51) 

(2.52) 

(2.53) 


m 

33  j 1  ^  3 


(2.54) 


Remark  2.8.  Since  the  distance  between  the  Chebysliev  collocation  points  xi, . . .  ,xm  near  the  ends  of 
the  interval  [—1, 1]  grow  as  to2  (see,  for  example,  [23]),  (2.47)  involves  division  by  small  numbers  when 
to  is  large.  As  discussed  in  [17],  this  is  remedied  by  using  the  trigonometric  identities 


—  Xj  =  2  sin  (  (i  +  j  —  2) 


2  (to  —  1) 


sin  (  ( j  -  i) 


2  (to.  -  1) 


>  *7 lj, 


(2.55) 


to  rewrite  (2.47)  into  the  form 


Di  i  =  — — 

hJ  2  a 


(— i  y+j 


3  sin  (  (i  +  j  —  2) 


2(m—  1) 


sm  ((*  -  j)  2(7^1)) 


i¥=j- 


(2.56) 


In  the  numerical  experiments  of  this  paper,  (2.56)  is  used  instead  of  (2.47)  to  compute  the  off-diagonal 
entries  of  D.  The  off-diagonal  entries  of  D  are  computed  using  a  similar  reformulation  of  (2.51)-(2.53). 
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3  Numerical  differentiation  via  PSWFs 


3.1  Differentiation  of  bandlimited  functions 

For  any  integer  k  >  1,  the  kth  order  derivative  of  a  bandlimited  function  /  :  [—  1, 1]  — >  R.  can  be  computed 
numerically  by  approximating  /  by  a  linear  combination  of  PSWFs,  and  then  differentiating  the  PSWFs 
k  times.  In  this  subsection,  we  construct  the  numerical  differentiation  scheme  for  computing  the  second 
derivative  of  /  using  the  PSWFs.  The  other  derivatives  of  /  can  be  computed  in  the  exact  same  manner. 

Let  /  :  [—1,1]  ->  R  be  a  bandlimited  function  with  bandlimit  c,  and  fix  some  small  e  >  0.  We  first 
consider  the  approximation  of  /  by  the  linear  combination  of  the  first  n  PSWFs: 

n—  1 

(3-1) 

3=0 

where  n  is  chosen  to  be  the  smallest  integer  such  that  the  corresponding  eigenvalue  of  the  operator  Fc 
has  magnitude  smaller  than  e,  i.e. , 

|A„_i|  >e>  |A„|.  (3.2) 

Next,  we  compute  the  coefficients  ctj  in  (3.1).  First,  we  apply  the  algorithm  of  [69]  to  construct  a 
quadrature  Xi, ,  xm,  w\, ... ,  wm  that  integrates  exactly  all  pairwise  products  of  ipo, ... ,  ipn-i  (see 
Remark  2.6).  Then,  we  evaluate  the  values  /i,  •  •  • ,  fm  of  /  at  the  nodes  Xi, . . . ,  xm,  and  apply  formula 
(2.36)  to  obtain 

a  =  P*WF,  (3.3) 

where 

a  =  (a0)  •  •  •  )  CCn- l)T, 

W  is  an  m  x  to  diagonal  matrix  with  entries  IV,;, t  =  y/wi,  and  P  is  an  m  x  n  matrix  with  entries 

Pi,j  =  \/w~i^j-i{xi),  i  =  1,  •  •  •  ,to;  j  =  1, . . .  ,n.  (3.6) 

Once  a  is  computed,  we  can  approximate  the  value  of  f"  on  any  x  £  [—1,1]  by  the  formula: 

n—  1 

(3J) 

3=0 

using  a  numerical  scheme  for  evaluating  tp"  at  arbitrary  points  (see  [69] ) . 

In  particular,  the  in  x  to  matrix  D  defined  by  the  formula 


(3.4) 

(3.5) 


D  =  UP*W. , 


(3.8) 


16 


where  U  is  an  m  x  n  matrix  with  entries 


Ui,j  =i>"-1(xi),  i  =  ljv. .  ,m;  j  =  1, . . .  ,n,  (3.9) 

is  a  linear  operator  taking  the  values  f\, ... ,  frn  of  the  function  /  at  the  interpolation  nodes  xi, . . . ,  xm  to 
the  values  di, . . . ,  dm,  where  dj  is  the  approximation  to  f"(xj)  via  (3.7).  We  refer  to  D  as  a  differentiation 
matrix  constructed  using  the  PSWFs.  In  Section  4.1.1,  we  demonstrate  the  accuracy  property  of  the 
above  scheme  via  several  numerical  experiments. 

Remark  3.1.  For  the  algorithm  of  [69],  the  nodes  of  the  quadrature  constructed  for  bandlimit  2c  and 
precision  eq  are  precisely  the  roots  of  ipf,  where  l  is  roughly  half  of  the  number  of  PSWFs  if?0  with 
eigenvalues  greater  than  eq.  Therefore,  the  quadrature  nodes  constructed  by  the  algorithm  of  [69],  and 
hence  the  interpolation  nodes  on  which  the  differentiation  matrix  D  is  constructed,  all  lie  in  the  interior 
of  the  interval  [—1, 1]. 

Remark  3.2.  Clearly,  the  scheme  described  above  is  not  restricted  to  functions  defined  on  [— 1, 1].  For 
example,  consider  a  function  g  :  [a,  6]  — >  R.  defined  by  transforming  another  function  /  :  [—  1, 1]  — >  R.  by 
the  formula 

g{x )  =  /(( 2x  -a-  b)/{b  -  a)),  x  £  [a,  b ];  (3.10) 

in  particular,  if  /  is  a  periodic  function  with  K  wavelengths  on  [—1,1],  then  g  is  also  periodic,  with 
K  wavelengths  on  [a,  b\.  Suppose  a  differentiation  matrix  D  (3.8)  that  approximates  /"  at  the  nodes 
Xi, ,  xm  £  [—1,1]  is  constructed  using  the  above  scheme,  then  the  matrix  D  obtained  by  the  formula 

Dij  =  (2/(6  —  a))2Dij,  i,j  =  1, . . .  ,m,  (3.11) 

approximates  g"  at  the  nodes  yi,.  ■■  ,ym  €=  [a,  b]  defined  by 

y%  =  (b-  a)xi/ 2  +  (a  +  b)/ 2,  i  =  1, . . . ,  m,  (3.12) 

to  the  same  order  of  relative  accuracies.  For  the  case  of  kth  derivative,  the  same  argument  applies  with 
the  exponent  2  in  (3.11)  replaced  by  k. 

3.2  Differentiation  matrices  incorporating  boundary  conditions 

For  almost  all  kinds  of  PDEs  encountered  in  practice,  boundary  conditions  are  needed  to  guarantee  the 
uniqueness  of  solutions.  The  numerical  solutions  of  these  PDEs  usually  require  the  incorporation  of  the 
boundary  conditions  into  the  discretized  systems.  One  approach  of  incorporating  boundary  conditions 
for  spectral  (and  pseudospectral)  methods  is  to  modify  the  interpolation  functions  to  make  them  satisfy 
the  boundary  conditions;  another  approach  is  to  add  additional  equations  to  the  discretized  systems 
to  enforce  the  boundary  conditions,  without  modifying  the  interpolation  functions  (see,  for  instance, 
[8,  23,  30,  63]).  In  this  subsection,  we  adopt  the  first  approach,  and  modify  the  scheme  in  Section  3.1  to 
construct  differentiation  matrices  that  incorporate  boundary  conditions. 
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As  a  specific  example,  we  consider  the  wave  equation  with  the  homogeneous  Dirichlet  boundary 
condition: 

Utt.  'Uxxi  X  ^  (  1)  1);  t  P>  0, 

u(x,0)  =  f(x),  x  £  [—1,1],  (3.13) 

u(— 1,  t)  =  u(l,t)  =  0,  t  >  0, 

and  construct  a  differentiation  matrix  D  that  discretizes  the  derivative  operator  uxx  and  incorporates 
the  boundary  condition.  The  main  objective  is  to  construct  using  the  PSWFs  an  orthonormal  set  of 
interpolation  functions  <pi, ...  ,4>k  in  L2([—  1, 1])  such  that  each  (f>j  satisfies  the  zero  boundary  condition: 

<A?(-1)  =  =  0,  j  =  l,...,k.  (3.14) 

First,  we  fix  a  bandlimit  c  and  precision  s,  and  take  the  first  n  PSWFs  . . . ,  V’n-i,  with  n  chosen  to  be 
the  smallest  integer  such  that  |An|  <  e.  We  then  subtract  from  each  ipj  a  linear  function  pij  interpolating 
the  values  of  tpj  at  the  endpoints  of  [—1,1]: 

H  0*0  =  ^(-i)  +  (V’j(i)  -  (3-15) 

The  resulting  set  of  functions 

Vb  =  Vb  “  MT  i  =  0,  (3.16) 

approximate,  with  an  error  of  the  order  of  e,  all  functions  /  :  [—1, 1]  — >  R.  that  have  bandlimit  c  and 
satisfy  the  zero  boundary  condition  /(— 1)  =  /( 1)  =  0. 

Next,  we  construct  an  orthonormal  basis  for  the  span  of  ipo, . . .  ,ipn- 1.  We  first  apply  the  algorithm 
of  [69]  to  construct  a  quadrature  Xi, . . . ,  xm,  w i, . .  • ,  wm  that  discretizes  ipo, . . . ,  ipn-i  (see  Remark  3.6 
below).  Then,  we  apply  the  pivoted  Gram-Schmidt  algorithm  (with  re-orthogonalization)  described  in 
[7]  to  the  in  x  n  matrix  A  with  entries 

Mxi)y/m  1pl(xi)s/Wl  ■■■  ^n-l{xi)y/wl 
1po{X2)y/W2  tpl(x2)^/W2  ■■■  ’lpn-l{X2)y/W2 

1^0  {Xm  )  y/  Xrn  Fl  (Xm  )  \fy'm  •••  l(.X  m)y/^rn 

obtaining  an  m  x  n  matrix  A.  Since  ipo, . . . ,  ipn~ \  approximates  all  bandlimited  functions  of  bandlimit  c 
to  precision  e,  and  i/j0 , . . . ,  V’n-i  are  obtained  by  subtracting  from  ipo, . . . ,  ipn_ i  two  degrees  of  freedom, 
the  matrix  A  has  £-rank  k  =  n  —  2.  We  take  the  first  k  columns  of  A  and  define  the  interpolation 
functions  <p\,  ■  ■  ■ ,  (j>k  by  their  values: 

(pj{xi)  =  Aij/y/vTi,  i  =  1,. . .  ,m,  j  =  1, . . . ,  k  (3.18) 

at  the  quadrature  nodes  xi, . . . ,  xm. 

Remark  3.3.  Although  <f>i, . . . ,  tpk  are  defined  by  their  values  on  the  quadrature  nodes  xi, . . . ,  xm,  their 
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values  on  any  arbitrary  point  x  €  [—1,1]  are  completely  determined  by  (3.16)  and  the  pivoted  Gram- 
Schmidt  procedure.  Specifically,  let  L  be  the  n  x  n  linear  transformation  corresponding  to  the  pivoted 
Gram-Schmidt  procedure: 

A  =  AL,  (3.19) 

then  the  value  of  is  obtained  by  applying  L  from  the  right  to  the  row  vector 

(3.20) 

and  taking  the  ith  element  of  the  resulting  row  vector.  In  particular,  each  of  the  0j  satisfies  the  boundary 
condition  (3.14). 

The  remaining  part  of  the  construction  is  similar  to  that  of  Section  3.1.  Let  /  :  [—1, 1]  — >  K  be  a 
function  of  bandlimit  c  such  that  /(— 1)  =  /( 1)  =  0,  then  /  can  be  approximated  by  a  linear  combination 
Of  0i,  ...  ,  (j>k ■ 

k 

'y  ]  ca j 4>j ,  (3.21) 

3=1 

with  an  error  of  the  order  of  e.  We  use  the  quadrature  nodes  Xi,... , xrn  as  interpolation  nodes,  and 
compute  the  m  x  m  differentiation  matrix  D  taking  the  values  (f(x i), . . . ,  /( xm))  to  the  approximating 
values  (di , . . . ,  dm  )  of  (/"(* i)  , . . . ,  f"(xm))  by  the  formula 

D  =  UP*W. ,  (3.22) 

where  W  is  an  m  x  m  diagonal  matrix  with  diagonal  entries  Wqy  =  yjwl,  P  is  an  m  x  k  matrix  with 
entries 

P,.,j  =  (t>j(xi)y/wi,  i  =  1, . . .  ,m;  j  =  1, . . .  ,k,  (3.23) 

and  U  is  an  m  x  k  matrix  with  entries 

Uij  =  (j)"{xi),  i  =  1, . . .  ,m;  j  =  1,  . . . ,  k.  (3.24) 

With  suitably  chosen  c  and  e,  the  differentiation  matrix  D  can  be  combined  with  a  suitable  time¬ 
marching  scheme  to  solve  (3.13)  to  any  desired  accuracy. 

Remark  3.4.  The  entries  (3.24)  of  U  are  obtained  by  applying  the  matrix  L  in  (3.19)  from  the  right  to 
the  m  x  n  matrix  S  defined  by 

Si,j  =  ipj-iixi),  i  =  1,  ••  •  ,m,  j  =  1, . . .  ,n,  (3.25) 

and  then  taking  the  first  k  columns  of  SL.  In  particular,  the  matrices  A  and  S  are  passed  together  to 
the  pivoted  Gram-Schmidt  algorithm  during  actual  implementation  of  the  scheme. 

Differentiation  matrices  that  incorporate  other  boundary  conditions,  such  as  the  Neumann  boundary 
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condition: 


ux(—l,t)  =  0,  u(l,t)=0,  t>  0, 


(3.26) 


are  similarly  constructed  by  modifying  the  functions  /ij  in  (3.15).  Likewise,  differentiation  matrices  for 
computing  first  derivatives  are  constructed  with  the  straightforward  modification  of  the  above  scheme. 
In  Section  4.2,  we  demonstrate  the  performance  of  the  scheme  by  applying  the  differentiation  matrices 
constructed  to  the  solution  of  several  time-dependent  PDEs  and  the  associated  eigenvalue  problems. 

Remark  3.5.  For  any  time  T  >  0,  we  combine  the  differentiation  matrix  D  with  a  time-marching  scheme 
to  solve  (3.13),  obtaining  numerical  approximation  (u\,...,um)  to  u(-,T)  at  the  nodes  xi, . . .  ,xm. 
Numerical  approximation  to  u(-,T )  at  another  set  of  nodes  yi,...,yi  are  then  computed  easily  by  an 
interpolation  scheme  based  on  (3.21)  and  (2.14). 

Remark  3.6.  In  order  to  construct  <f>i, . . .  ,(f>k  by  the  Gram-Schmidt  algorithm,  the  quadrature  xi, . . . ,  xm, 
Wi . . . ,  wm  has  to  be  chosen  such  that  it  discretizes  ipo, . . . ,  V’n-i-  In  the  numerical  experiments  in 
Section  4.1.2,  we  construct  using  the  algorithm  of  [69]  a  quadrature  integrating  all  functions  on  [—1, 1]  of 
bandlimit  2c  to  precision  e2.  Other  quadratures,  such  as  the  Gaussian-Legendre  quadrature,  can  also  be 
used.  As  observed  in  [69],  however,  the  latter  choice  requires  more  nodes  to  achieve  the  same  accuracy. 

Remark  3.7.  The  non-zero  eigenvalues  of  the  differentiation  matrix  D  constructed  by  the  above  scheme 
depends  only  on  the  PSWFs  tpo, ... ,  ipn-h  and  hence  only  on  the  choices  of  c  and  e.  In  particular,  they 
do  not  depend  on  the  particular  quadrature  used  or  its  number  of  nodes  m,  as  long  as  it  is  accurate 
enough  to  discretize  the  functions  ipo, . . . ,  V’n-i-  In  other  words,  given  fixed  c  and  e  and  an  appropriately 
chosen  quadrature,  the  non-zero  part  of  the  spectrum  of  D  is  independent  of  its  dimension  m. 

4  Numerical  results 

4.1  Accuracy  and  stability  properties 

In  this  subsection,  we  demonstrate  results  regarding  the  accuracy  and  stability  properties  of  the  numerical 
differentiation  schemes  described  in  Section  3.  The  schemes  were  implemented  in  FORTRAN  77.  All 
implementations  were  carried  out  in  double  (16-digit)  precision  arithmetic  except  otherwise  noted. 

4.1.1  Accuracies  of  the  differentiation  matrices 

In  this  subsection,  we  present  results  related  to  the  accuracy  of  the  differentiation  matrices  constructed 
by  the  algorithm  described  in  Section  3.1.  Let  D  be  the  m  x  m  differentiation  matrix  constructed 
using  the  algorithm  with  bandlimit  c  and  precision  e,  taking  the  values  of  functions  /  :  [—1,1]  — >  R 
at  the  interpolation  nodes  xi,...,xm  €  [—1,1]  to  its  approximate  derivatives  at  the  nodes.  For  each 
a  >  0,  we  apply  D  to  obtain  approximations  of  the  derivatives  of  the  functions  sin(acc)  and  cos  (ax)  at 
x\, . . .  xm.  We  denote  the  approximations  by  Ui,  u^, . . . ,  «2m  and  the  corresponding  analytical  derivatives 
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by  Vi,  V2,  ■  ■  ■ ,  V2 mi  and  measure  the  errors  using  the  following  three  quantities: 


Ea,i  — 

max  I Ui  —  vA. 

(4.1) 

l<i<2m 

Ea,  2  = 

V  2  TO 

(4.2) 

Ea,  3  = 

/eS iK-«i)2 

V  E&ttf  • 

(4.3) 

In  other  words,  Ea_\,  Ea^ 2,  and  Ea ^  are  the  maximum  absolute  errors,  root-mean-squared  errors,  and 
relative  I2  errors  respectively.  Table  2  shows,  for  matrices  D  constructed  using  e  =  10“ 14  and  varying 
values  of  c,  the  maximum  errors  maxog[0  c]  Ea  i  and  maxa6[0c]  Ea  2  over  a  large  number  of  values  of  a 
sampled  on  [0 ,  c] ,  as  well  as  the  relative  I2  error  ECt 3.  Table  3,  on  the  other  hand,  shows  the  above 
errors  for  a  fixed  bandlimit  c  =  1000  and  varying  values  of  e.  In  the  tables,  m  denotes  the  number 
of  nodes  x\, . . .  ,xm  used  in  the  construction  of  D ,  and  hence  the  dimension  of  D,  while  n  denotes  the 
number  of  interpolating  PSWFs  -i/’j  corresponding  to  the  choice  of  c  and  £  (see  (3.1)-(3.2)).  Figure  1 
shows  the  absolute  errors  across  the  nodes  on  [—1, 1]  when  D  is  used  to  differentiate  the  single  function 
f(x)  =  sin(ca:),  where  D  is  constructed  with  £  =  10~14  and  c  =  50,1200.  Lastly,  Figure  2  shows,  for 
D  constructed  using  c  =  50,100  and  £  =  10-8,10-14,  the  errors  Ea^\  across  a  large  number  of  values 
of  a  sampled  on  [0,200].  Results  associated  with  both  first  and  second  derivatives  are  shown  in  the 
above  tables  and  figures.  In  the  following  discussion,  we  further  explain  our  results  and  make  several 
observations. 

1.  In  these  experiments,  the  nodes  aq, . . . ,  xm  are  chosen  to  be  the  nodes  of  a  quadrature  integrating 
functions  of  bandlimit  2c  to  precision  e2  x  10~ 10 ,  constructed  using  the  algorithm  of  [69].  This 
ensures  that  the  columns  of  the  matrix  P  (3.6)  are  orthonormal,  to  machine  precision.  From 
inspection  of  Theorem  2.2,  the  dependence  of  the  number  of  PSWFs  ipj  on  the  eigenvalue  cutoff 
£  is  weak.  Hence,  choosing  the  precision  £2  x  10~10  leads  to  only  slightly  larger  numbers  of 
interpolation  nodes  compared  to  those  used  in  the  numerical  experiments  of  [69]  (see  Remark  3.1). 

2.  As  expected  again  from  Theorem  2.2,  the  dependence  of  the  number  of  nodes  m  on  the  precision 
£  is  weak.  For  example,  Table  3  shows  that  for  c  =  1000,  we  roughly  gain  one  digit  of  accuracy  by 
increasing  the  number  of  nodes  m  by  about  3. 

3.  From  Table  2,  we  see  that  the  sampling  factor ,  defined  as  the  ratio  of  m  over  n,  is  close  to  1,  and 
approaches  closer  to  1  as  c  increases.  The  reason  is  because  m  is  roughly  half  of  the  number  of 
PSWFs  ip?0  with  eigenvalues  greater  than  £2  x  10-10  (see  Remark  3.1),  and  the  latter  quantity  is 
roughly  double  that  of  n,  due  to  Theorem  2.2. 

4.  For  large  c,  the  ratio  of  the  number  of  nodes  m  over  ^  is  close  to  two,  regardless  of  the  accuracy 
requirement.  This  means  that,  for  large  c,  a  slightly  more  than  two  points  per  wavelength  are 
needed  to  achieve  any  desired  accuracy. 
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5.  For  each  a  >  0,  the  maximum  absolute  error  Eai  always  occurs  near  the  ends  of  the  interval  [— 1, 1], 
regardless  of  the  differentiation  matrix  used  to  approximate  the  derivatives.  Figure  1  shows,  as  an 
example,  that  for  first  derivatives  the  accuracies  are  usually  1  —  2  digits  lower  near  the  endpoints 
compared  to  the  rest  of  the  interval,  while  for  second  derivatives  the  accuracies  are  usually  3  —  4 
digits  lower  near  the  endpoints.  This  is  attributed  to  the  fact  that  the  magnitudes  of  the  derivatives 
of  each  PSWF  ijjj(x)  rises  sharply  as  x  approaches  the  ends  of  [—1,1],  which,  in  turn,  implies  that 
the  norms  of  the  first  and  last  few  rows  of  the  differentiation  matrices  are  usually  several  orders 
of  magnitude  larger  than  those  of  the  rows  near  the  middle.  As  an  example,  Table  4  shows  the  I2 
norms  of  the  first  and  middle  rows  of  the  differentiation  matrices  constructed  using  different  values 
of  c  and  e. 

6.  From  Table  2,  we  see  that  the  absolute  measures  of  errors  maxag[o,ol  Eat\  and  maxaeroiC]  Ea<2 
increase  roughly  at  a  rate  of  c  and  c2  for  first  and  second  derivatives,  respectively.  This  is  consistent 
with  the  fact  that  the  amplification  factors  for  taking  the  first  and  second  derivatives  of  the 
functions  sinfca:)  and  cos(c;r)  are  about  c  and  c2,  respectively.  For  example,  we  see  from  Figure  1 
that  for  c  =  1200  and  £  =  10-14,  we  get  an  absolute  error  of  about  10”8  near  the  center  of  [—1, 1] 
when  approximating  the  second  derivatives  of  f(x)  =  sin(1200a’),  losing  about  6  digits  relative  to 
the  prescribed  e. 

7.  From  Figure  2,  we  see  that  the  maximum  absolute  error  Ea  \  for  a  differentiation  matrix  constructed 
using  bandlimit  c  is  almost  uniform  over  all  0  <  a  <  c.  The  accuracy  decreases  exponentially  once 
a  increases  beyond  the  prescribed  bandlimit  c  of  the  matrix,  with  the  rate  of  decrease  almost 
the  same  for  c  =  50  and  c  =  100.  Our  more  extensive  tests  show  that  for  £  =  10” 14  and  any 
prescribed  c,  no  accuracy  is  obtained  when  a  gets  about  5  beyond  c.  This  amounts  to  fewer  than 
two  wavelengths  over  the  interval  [—1, 1]. 

8.  From  Tables  2  and  3  we  see  that,  given  any  c,  the  relative  I2  errors  ECt 3  are  consistently  about 
1  and  2  digits  less  than  the  prescribed  e,  for  first  and  second  derivatives  respectively.  Figure  3 
shows  the  errors  Eai 3  across  0  <  a  <  c  for  the  second  derivative  matrix  constructed  with  c  =  1000 
and  £  =  10”8, 10-14,  indicating  that  Ea  3  increases  as  a  decreases  towards  0.  This  is  expected, 
given  the  observation  in  (7)  that  the  matrices  constructed  by  the  algorithm  give  roughly  uniform 
absolute  errors  for  all  0  <  a  <  c.  Therefore,  it  is  optimal  in  terms  of  both  absolute  and  relative 
errors  to  approximate  the  derivatives  of  functions  /  :  [—1,1]  — >  R.  with  bandlimit  roughly  c  with 
differentiation  matrices  constructed  using  c. 
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Table  4:  I2  norms  of  the  rows  of  the  differentiation  matrices  constructed  for  first  and  second  derivatives 
using  varying  values  of  c  and  e.  E\  denotes  the  I2  norm  of  the  first  row  of  the  matrix,  while  E2  denotes 
the  I2  norm  of  the  jth  row,  where  j  is  the  largest  integer  less  than  or  equal  to  in/2. 


c 

£ 

m 

j  First  derivatives  I 

Second  derivatives 

E 1 

e2 

Ex 

E2 

50 

1.0E-08 

57 

6.84E+02 

3.42E+01 

2.32E+05 

1.60E+03 

50 

1.0E-14 

66 

1.26E+03 

3.84E+01 

8.11E+05 

2.12E+03 

200 

1.0E-08 

161 

2.59E+03 

1.20E+02 

3.25E+06 

1.96E+04 

200 

1.0E-14 

173 

4.77E+03 

1.26E+02 

1.14E+07 

2.14E+05 

600 

1.0E-08 

422 

7.50E+03 

3.51E+02 

2.69E+07 

1.66E+05 

600 

1.0E-14 

437 

1.35E+04 

3.56E+02 

8.92E+07 

1.71E+05 

1400 

1.0E-08 

936 

1.85E+04 

8.13E+02 

1.66E+08 

8.90E+05 

1400 

1.0E-14 

954 

3.23E+04 

8.18E+02 

5.15E+08 

9.07E+05 

In  addition  to  looking  at  the  accuracies  of  the  differentiation  matrices  on  functions  of  the  form  sin(aa;) 
and  cos(aa;),  we  look  at  the  accuracies  when  the  matrices  are  used  to  approximate  the  second  derivatives 
of  the  functions  exp(— 1000a;2),  sin(207ra;  +  sin(207rx)),  and  cos2(300a;)  on  the  interval  [—1, 1].  Figure  4 
shows  the  relative  I2  errors,  with  the  differentiation  matrices  constructed  using  e  =  10_  14  and  varying 
values  of  c  in  the  interval  [0, 1000].  The  errors  are  computed  on  the  nodes  corresponding  to  the  matrices. 
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Figure  2:  The  errors  Ea< \  across  a  £  [0,  200]  for  (a)  first  and  (b)  second  derivatives,  where  D  is  constructed 
using  c  =  50, 100  and  e  =  10-8, 10”14. 
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a 


Figure  3:  The  errors  E0t 3  across  a  £  [0, 1000]  for  the  second  derivative  matrix  constructed  using  c  =  1000 
and  e  =  10-8, 10~14. 


f(x)  =  e~1000x2  — ^ 
f(x)  =  sin(207ra;  +  sin(207ra’))  — *- 
f(x)  =  cos2  (300a:)  — & 


c 


Figure  4:  Relative  I2  errors  of  second  derivatives  for  the  functions  /( x)  =  exp(— 1000a:2),  sin(207ra:  + 
sin(207ra;)),  and  cos2(300a;)  defined  on  [—1,1].  The  second  derivative  matrices  are  constructed  using 
£  =  10“ 14  and  varying  values  of  c  in  [0, 1000]. 
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4.1.2  Eigenvalues  of  the  differentiation  matrices 

One  of  the  main  applications  of  differentiation  matrices  is  the  solution  of  time-dependent  PDEs.  Often, 
the  spatial  derivatives  of  a  time-dependent  PDE  are  discretized  by  differentiation  matrices  that  incor¬ 
porate  the  boundary  conditions;  the  resulting  semi-discretized  system  is  then  solved  using  a  numerical 
ODE  solver,  such  as  the  Runge-Kutta  scheme  with  a  suitably  chosen  time-step  At.  In  many  practical 
situations,  the  stability  of  such  a  scheme  for  the  solution  of  time-dependent  PDEs  is  determined  by  the 
eigenvalues  of  the  discretized  spatial  operator.  Namely,  the  eigenvalues  of  the  discretized  spatial  oper¬ 
ator,  multiplied  by  the  time-step  At,  have  to  lie  inside  the  stability  region  of  the  ODE  solver  in  order 
to  ensure  stability.  In  this  subsection,  we  present  results  related  to  the  eigenvalues  of  the  differentiation 
matrices  constructed  using  the  scheme  described  in  Section  3.2,  and  compare  them  to  the  eigenvalues 
of  the  differentiation  matrices  constructed  using  the  Cliebyshev  collocation  method.  Unless  otherwise 
stated,  all  eigenvalues  are  computed  numerically  in  double  precision  arithmetic. 

Second-order  differentiation  matrices 

We  first  consider  the  wave  equation  with  the  homogeneous  Diriclilet  boundary  condition: 

Utt.  'U'xxt  X  ^  (  1?  1);  t  0, 

u(x,0)  =  f(x),  x  £  [—1, 1],  (4.4) 

u(—  l,t)  =  =  0,  t  >  0, 

and  construct  using  the  scheme  described  in  Section  3.2  a  differentiation  matrix  D  that  discretizes  the 
operator  uxx  and  incorporates  the  boundary  condition.  The  constructed  matrix  D  also  discretizes  the 
eigenvalue  problem 


Uxx  —  A'U,  X  (E  (  1,1), 
u(±l)  =  0, 


(4.5) 


which  has  analytical  eigenvalues 

j 27T2 

Aj  =  — ,  j  =  1,2, , 

and  corresponding  eigenfunctions 


(4.6) 


u(x)  =  sin  J7r(^+  j  =  l,  2, ... .  (4.7) 

Table  5  shows  the  non-zero  eigenvalues  of  the  matrix  D  constructed  with  c  =  20-7T  and  e  =  10-13.  There 
are  66  non-zero  eigenvalues,  which  is  the  same  as  the  number  of  functions  k  in  the  orthonormal  set 
<f>i, . . . ,  ff>k  (3.18)  corresponding  to  the  chosen  c  and  e.  For  ease  of  comparison  with  (4.6),  the  eigenvalues 
Ai, . . . ,  Xqq  in  Table  5  are  scaled  by  the  factor  4/7T2.  From  the  table,  we  see  that  Ai, . . . ,  A66  are  real, 
distinct,  and  negative;  and  Ai, . . . ,  A40  are  accurate  to  at  least  13  digits  with  respect  to  the  eigenvalues 
of  the  continuous  problems  (4.6),  which  is  expected  given  the  choice  of  the  c  used  in  the  construction  of 


29 


Figure  5:  Plot  of  the  absolute  values  of  the  eigenvalues  Ai, . . . ,  A66  in  Table  5  compared  with  the  analytical 
eigenvalues  (4.6)  drawn  on  the  solid  line. 

D.  For  k  >  40,  the  magnitudes  of  A j  start  to  grow  exponentially,  and  they  no  longer  approximate  the 
eigenvalues  (4.6).  Figure  5  shows  the  plot  of  | A;  compared  with  the  eigenvalues  (4.6). 

In  addition,  we  compute  the  eigenvalues  of  the  matrix  D  constructed  using  c  =  6007T  and  e  =  10~13. 
There  are  1246  non-zero  eigenvalues  A  j.  All  of  them  are  again  real,  distinct,  and  negative,  and  agree 
with  (4.6)  to  at  least  13  digits,  for  all  1  <  j  <  1200.  Figure  6  shows  the  plot  of  | Ay  |  compared  with  (4.6). 
The  plot  shows  that  |  A  3  |  begin  to  grow  exponentially  for  j  >  1200. 

Remark  4.1.  We  construct  the  above  differentiation  matrices  using  the  scheme  described  in  Section  3.2, 
in  which  we  use  the  algorithm  of  [69]  to  construct  a  quadrature  integrating  all  functions  on  [—1,1]  of 
bandlimit  2c  to  precision  e2 .  Our  numerical  experiments  indicate  that  such  choice  of  quadrature  is 
sufficient  in  discretizing  the  functions  i/j o, . . .  ,ipn-i  (3.16).  Unless  otherwise  stated,  we  will  adopt  this 
choice  of  quadrature  for  the  rest  of  this  paper.  As  pointed  out  in  Remark  3.7,  as  long  as  the  choice  of  c 
and  e  is  fixed,  increasing  the  precision  of  the  quadrature  beyond  £  ,  and  hence  increasing  the  size  of  D1 
has  no  effect  on  the  non-zero  eigenvalues. 

In  the  numerical  solution  of  a  time-dependent  PDE  using  a  combination  of  a  differentiation  matrix 
with  an  explicit  time-marching  scheme,  the  largest  size  of  the  time-step  At  that  maintains  stability  is 
often  determined  by  the  spectral  radius  of  the  differentiation  matrix,  defined  as  the  maximum  moduli 
of  its  eigenvalues.  In  the  following,  we  look  at  the  spectral  radii  of  the  differentiation  matrices  D  dis¬ 
cretizing  (4.5)  constructed  using  the  scheme  described  in  Section  3.2,  and  compare  them  with  those  of 
the  differentiation  matrices  constructed  using  the  Chebyshev  collocation  method.  Let  Dc  &  lSlNxN  de¬ 
note  the  second  derivative  matrix  constructed  using  the  Chebyshev  collocation  method  on  N  collocation 
points,  as  described  in  Section  2.3.2.  There  are  several  ways  of  imposing  the  zero  boundary  conditions 
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Table  5:  Non-zero  eigenvalues  Ai, . . . ,  Ag6  of  the  differentiation  matrix  D  discretizing  the  problem  (4.5), 
constructed  by  the  scheme  described  in  Section  3.2  using  c  =  207t  and  e  =  10-13.  Each  A j  is  scaled  by 
the  factor  4/-7T2. 


-.1000000000000189E+01 
-.9000000000000133E+01 
-.2500000000000034E+02 
-.4900000000000004E+02 
-.8100000000000016E+02 
-.1210000000000001E+03 
-.1690000000000004E+03 
-.2250000000000000E+03 
-.2889999999999992E+03 
-.3610000000000005E+03 
-.4409999999999986E+03 
-.5290000000000006E+03 
-.6250000000000019E+03 
-.7289999999999994E+03 
-.8409999999999986E+03 
-.9609999999999964E+03 
-.  1089000000000000E+04 
-.1224999999999996E+04 
-.1369000000000002E+04 
-.1520999999999995E+04 
-.1681000000000003E+04 
-.  1849000000000964E+04 
-.2025000012389790E+04 
-.2209009537738922E+04 
-.2401989463349080E+04 
-.2620587546324507E+04 
-.2930504310944726E+04 
-.3435568410972565E+04 
-.4308191253136334E+04 
-.5951192152853712E+04 
-.9541230907887420E+04 
-.1984784397165518E+05 
-.7558406595445475E+05 


-.4000000000000032E+01 
-.1600000000000005E+02 
-.3600000000000038E+02 
-.6400000000000021E+02 
-.1000000000000002E+03 
-.1440000000000000E+03 
-.1960000000000000E+03 
-.2559999999999995E+03 
-.3240000000000014E+03 
-.4000000000000004E+03 
-.4840000000000009E+03 
-.5760000000000006E+03 
-.6759999999999987E+03 
-.7840000000000028E+03 
-.9000000000000000E+03 
-.1024000000000003E+04 
- .  1 156000000000000E+04 
-.1296000000000005E+04 
-.1444000000000000E+04 
-.1600000000000006E+04 
-.1764000000000001E+04 
-.1936000000030938E+04 
-.  21 16000113341 137E+04 
-.2304041387127145E+04 
-.  25024905471 17806E+04 
-.2737428098562087E+04 
-.3083256186436362E+04 
-.3645848773066469E+04 
-.4612338137669602E+04 
-.6426158975995367E+04 
-.1038294717510027E+05 
-.2173440910500743E+05 
-.8310762526925356E+05 
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Figure  6:  Plot  of  the  absolute  values  of  the  non-zero  eigenvalues  A j  of  the  differentiation  matrix  D 
discretizing  the  problem  (4.5),  constructed  by  the  scheme  described  in  Section  3.2  using  c  =  6007T  and 
e  =  10-13.  There  are  1246  of  these  eigenvalues,  and  they  are  compared  with  the  analytical  eigenvalues 
(4.6)  drawn  on  the  solid  line. 


u(±l)  =  0  on  Dc  (see,  for  instance,  [8,  23,  63]),  one  of  which  is  to  strip  Dc  of  its  first  and  last  rows  and 
columns,  resulting  in  the  (TV  —  2)  x  (TV  —  2)  matrix  Dc-  In  [29],  it  is  proven  that  all  eigenvalues  of  Dc 
are  real,  distinct,  and  negative. 

Tables  6  and  7  show  the  spectral  radii  of  D  constructed  with  varying  values  of  c  and  e  =  10~7 
and  10~13,  respectively.  All  computed  eigenvalues  are  either  zero,  or  real,  distinct,  and  negative,  and 
the  spectral  radii  are  compared  with  those  of  Dc-  The  following  notation  is  used  when  presenting  the 
results. 


to  =  number  of  nodes  used  in  the  construction  of  D,  with 
bandlimit  parameter  chosen  as  c 
N  =  number  of  nodes  used  in  the  construction  of  Dc, 
corresponding  to  Dc  of  dimension  N  —  2 
Ec  =  relative  1%  error  when  D  or  Dc  is  used  to  approximate 
the  second  derivatives  of  the  functions  sin(ca’)  —  x  sin(c) 
and  cos (cx)  —  cos(c)  on  [—1,1] 
p  =  spectral  radius  of  D  or  Dc 

The  errors  Ec  are  computed  on  the  quadrature  nodes  or  the  Chebyshev  collocation  points  corresponding 
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to  D  or  Dq  respectively;  and  the  functions  sin(ca’)  —  ccsin(c)  and  cos (cx)  —  cos(c)  are  just  sine  and 
cosine  functions  with  linear  functions  subtracted  to  make  them  satisfy  the  zero  condition  at  x  =  ±1. 
For  each  value  of  c,  we  choose  the  number  of  nodes  N  for  the  Chebyshev  differentiation  matrix  Dq  to 
give  similar  error  Ec  as  does  the  corresponding  D.  The  number  of  nodes  to  for  D  is  determined  by  the 
choice  of  c  and  e,  as  well  as  the  choice  of  quadrature  (see  Remark  4.1).  In  the  following,  we  make  several 
observations  and  remarks  regarding  the  results. 

1.  From  Tables  6  and  7,  we  observe  that  for  fixed  e,  the  spectral  radius  p  of  D  grows  like  c2.  Figure  7 
shows  the  plot  of  p  against  c  on  a  log-log  scale,  which  verifies  the  observation. 

2.  Since,  as  pointed  out  in  Section  4.1.1,  the  ratio  of  number  of  nodes  to  to  c/7r  is  roughly  equal  to 
2  for  large  c,  we  can  interpret  the  observation  of  (1)  as  saying  that  the  spectral  radius  of  D  grows 
roughly  as  m2  for  large  to.  However,  caution  must  be  taken  when  interpreting  the  results  this  way, 
since  to  is  not  an  independent  variable  in  our  construction  of  the  matrices. 

3.  The  spectral  radius  p  of  Dq,  on  the  other  hand,  grows  as  N 4,  as  illustrated  in  Figure  8.  This  is 
a  well-known  property  of  collocation  methods  on  both  Chebyshev  and  Legendre  points  (see,  for 
instance,  [68]),  and  it  imposes  strict  stability  condition  on  the  time-step  allowed  in  an  explicit 
time-stepping  scheme. 

For  example,  from  Table  7,  we  see  that  the  differentiation  matrix  D  constructed  with  c  =  2000 
and  e  =  10-13  has  spectral  radius  1.53E+08,  while  the  Chebyshev  matrix  Dq  that  yields  similar 
error  Ec  has  spectral  radius  9.19E+11.  This  means  that  for  the  solution  of  the  parabolic  equation 
ut  =  uxx,  the  step-size  required  to  maintain  stability  when  using  Dq  to  discretize  uxx  are  about 
6000  times  smaller  than  those  needed  when  using  D  to  discretize  uxx.  For  the  solution  of  the 
hyperbolic  equation  utt  =  uxx ,  the  ratio  of  the  required  step-sizes  is  about  V6000  «  77.  The 
numerical  results  in  Section  4.2  further  illustrate  the  difference  in  stability  requirements  implied 
by  the  spectral  radii  corresponding  to  the  two  differentiation  schemes. 

4.  From  Tables  6  and  7,  we  see  that  for  functions  with  large  bandlimits,  the  number  of  nodes  required 
by  the  Chebyshev  collocation  method  is  about  ir/2  times  more  than  the  nodes  required  by  the 
scheme  described  in  Section  3.2.  The  results  are  in  agreement  with  Theorem  2  of  [68],  which  asserts 
that  asymptotically  as  the  bandwidth  a  increases,  about  ir  points  per  wavelength  are  needed  for 
the  interpolation  of  the  functions  sin(a:r)  and  cos(a;r)  on  the  Chebyshev  collocation  points. 


33 


Table  6:  Spectral  radii  p  of  the  differentiation  matrices  D  discretizing  the  problem  (4.5),  constructed 
by  the  scheme  described  in  Section  3.2,  with  £  =  10  '  and  varying  values  of  c.  For  each  c,  the  spectral 
radius  is  compared  with  that  of  the  Chebyshev  matrix  Dc  that  gives  comparable  error  Ec. 


c 

m 

PSWFs 

Ec 

P 

N 

Chebyshev 

Ec  p 

5 

13 

2.25E-05 

5.63E+02 

15 

8.56E-06 

1.87E+03 

10 

18 

1.94E-06 

2.29E+03 

23 

2.13E-06 

1.12E+04 

20 

26 

2.61E-06 

6.93E+03 

36 

1.66E-06 

7.14E+04 

40 

41 

6.20E-06 

2.13E+04 

58 

7.16E-06 

5.01E+05 

50 

48 

1.95E-06 

3.77E+04 

71 

1.65E-06 

1.14E+06 

80 

68 

3.48E-06 

8.47E+04 

104 

2.02E-06 

5.33E+06 

100 

82 

1.94E-06 

1.40E+05 

126 

1.68E-06 

1.16E+07 

200 

150 

1.72E-06 

5.42E+05 

233 

1.22E-06 

1.37E+08 

400 

280 

1.21E-06 

2.19E+06 

441 

1.42E-06 

1.78E+09 

500 

344 

1.62E-06 

3.28E+06 

545 

9.81E-07 

4.15E+09 

800 

534 

2.40E-06 

7.94E+06 

850 

2.44E-06 

2.46E+10 

1000 

665 

2.34E-06 

1.23E+07 

1055 

1.67E-06 

5.85E+10 

1200 

793 

2.18E-06 

1.77E+07 

1258 

1.92E-06 

1.18E+11 

1600 

1049 

1.65E-06 

3.25E+07 

1664 

1.84E-06 

3.62E+11 

2000 

1302 

1.19E-06 

5.30E+07 

2070 

1.40E-06 

8.68E+11 

Table  7:  Spectral  radii  p  of  the  differentiation  matrices  D  discretizing  the  problem  (4.5),  constructed  by 
the  scheme  described  in  Section  3.2,  with  £  =  10-13  and  varying  values  of  c. 


c 

m 

PSWFs 

Ec 

P 

N 

Chebyshe 

Ec 

V 

p 

5 

18 

3.77E-11 

2.82E+03 

23 

2.33E-12 

1.12E+04 

10 

24 

1.942-11 

7.84E+03 

32 

1.33E-12 

4.40E+04 

20 

34 

6.26E-12 

2.55E+04 

46 

4.01E-12 

1.95E+05 

40 

50 

8.39E-12 

8.18E+04 

71 

1.21E-11 

1.14E+06 

50 

59 

4.53E-12 

1.28E+05 

84 

5.49E-12 

2.25E+06 

80 

81 

3.26E-12 

3.12E+05 

120 

3.46E-12 

9.50E+06 

100 

95 

3.09E-12 

4.74E+05 

143 

3.09E-12 

1.93E+07 

200 

161 

1.98E-11 

1.57E+06 

250 

1.87E-11 

1.82E+08 

400 

292 

1.12E-11 

6.38E+06 

467 

1.65E-11 

2.23E+09 

500 

360 

4.37E-12 

1.05E+07 

575 

3.70E-11 

5.14E+09 

800 

554 

3.68E-12 

2.68E+07 

890 

4.11E-11 

2.96E+10 

1000 

682 

6.61E-12 

4.02E+07 

1080 

2.38E-10 

6.42E+10 

1200 

811 

1.00E-11 

5.64E+07 

1290 

1.13E-10 

1.31E+11 

1600 

1067 

1.24E-11 

9.81E+07 

1690 

5.67E-10 

3.85E+11 

2000 

1324 

1.05E-11 

1.53E+08 

2100 

4.37E-10 

9.19E+11 
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First-order  differentiation  matrices 


In  the  remaining  part  of  this  subsection,  we  look  at  the  eigenvalues  of  the  first  derivative  matrices  con¬ 
structed  by  the  scheme  described  in  Section  3.2,  and  compare  them  to  those  of  the  matrices  constructed 
using  the  Chebysliev  collocation  method.  We  consider  the  first-order  hyperbolic  initial  boundary  value 
problem 

ut  =ux,  x  €  (-1, 1),  t  >  0, 

u(x,0)=f{x),  x  £  [—1, 1],  (4.8) 

u(l,t)  =  0,  t  >  0, 

and  discretize  the  spatial  derivative  operator  ux  by  the  obvious  modification  of  the  scheme  described  in 
Section  3.2.  Namely,  we  change  the  /q  defined  in  (3.15)  into  the  constant  functions 


=  i>i(  1),  (4.9) 

and  change  the  second  derivatives  in  (3.24)  and  (3.25)  into  first  derivatives.  We  denote  the  resulting 
differentiation  matrix  by  D. 

Remark  4.2.  The  other  modification  is  that  here  we  are  subtracting  from  the  PSWFS  ipo, . . . ,  ipn-i  one 
degree  of  freedom,  so  the  matrix  A  defined  in  (3.17)  has  e-rank  equal  to  n  —  1,  where  £  is  the  precision 
with  respect  to  which  ipi, . . .  ,%l>n  are  chosen. 

Figure  9  shows  on  the  complex  plane  the  eigenvalues  of  the  matrix  D  constructed  using  c  =  207t  and 
£  =  10^13.  There  are  67  non-zero  eigenvalues,  which  equals  the  number  of  functions  k  in  the  orthonormal 
set  (j>i, ...  ,(j>k  (3.18)  corresponding  to  c  and  £.  All  eigenvalues  of  D  lie  on  the  left  half-plane,  so  stability 
is  ensured  when  we  discretize  ux  in  (4.8)  by  D  and  solve  the  resulting  system  by  an  explicit  time-stepping 
scheme  with  sufficiently  small  time-step.  Most  of  the  eigenvalues  are  distributed  along  a  bow-shaped 
arc  around  the  origin,  with  a  few  outliers  extending  beyond  the  arc.  The  farthest  outlying  eigenvalues 
A±  =  (—0.407,  ±226)  almost  touch  the  imaginary  axis,  so  the  spectral  radius  of  D  is  roughly  equal  to 
the  magnitude  of  their  imaginary  part.  Table  8  shows  the  spectral  radii  of  D  constructed  using  varying 
values  of  c  and  £  =  10— 13 ,  and  compares  them  with  the  spectral  radii  of  the  matrices  Dc  constructed 
using  the  Chebyshev  collocation  method.  More  specifically,  Dc  is  obtained  by  constructing  the  first 
derivative  matrix  Dc  as  described  in  Section  2.3.2,  and  then  removing  its  first  row  and  column.  The 
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following  notation  is  used  when  presenting  the  results. 


in  =  number  of  nodes  used  in  the  construction  of  D1  with 
bandlimit  parameter  chosen  as  c 
N  =  number  of  nodes  used  in  the  construction  of  Dc, 
corresponding  to  Dq  of  dimension  N  —  1 
Ec  =  relative  1 2  error  when  D  or  Dq  is  used  to  approximate 
the  first  derivatives  of  the  functions  sin(cx)  —  sin(c) 
and  cos(ca’)  —  cos (c)  on  [—1, 1] 
p  =  spectral  radius  of  D  or  Dq 

For  each  value  of  c,  we  choose  the  number  of  nodes  N  for  Dc  to  give  similar  error  Ec  as  does  the 
corresponding  D.  In  addition,  we  plot  in  Figures  10  and  11  the  eigenvalues  of  D  and  Dc  corresponding 
to  several  selected  results  in  Table  8.  In  the  following,  we  make  several  observations  and  remarks  based 
on  the  presented  results,  and  on  the  results  of  our  more  extensive  experiments. 

1.  In  contrast  to  the  second-order  problem  (4.4),  the  first  derivative  operator  ux  in  (4.8),  with  the 
boundary  condition  u(l)  =  0,  does  not  correspond  to  any  eigenvalue  problem.  Therefore,  the 
eigenvalues  of  D  and  Dc  are  not  approximations  to  any  physically  meaningful  eigenvalues,  and 
are  entirely  numerical  in  origin. 

2.  All  eigenvalues  of  D  and  Dc  lie  on  the  left  half-plane.  Most  of  the  eigenvalues  of  D  are  distributed 
along  a  bow-shaped  arc,  with  some  outlying  eigenvalues  extending  along  the  imaginary  axis  (see 
Figure  10).  As  we  fix  e  and  increase  c,  the  imaginary  parts  of  the  outlying  eigenvalues  increase 
proportionately,  while  the  real  parts  of  all  eigenvalues  remain  bounded.  In  particular,  all  of  the 
eigenvalues  of  D  corresponding  to  the  results  in  Table  8  have  real  parts  bounded  by  20.  On  the 
other  hand,  as  observed  in  [64],  most  of  the  eigenvalues  of  Dc  are  distributed  along  an  arc  or  a 
loop,  again  with  some  outliers  (see  Figure  11).  In  contrast  to  D ,  both  the  real  and  imaginary  parts 
of  the  outlying  eigenvalues  grow  as  we  increase  N . 

3.  For  fixed  e.  the  spectral  radius  of  D  grows  as  c.  Again,  this  can  be  interpreted  as  that  the  spectral 
radius  of  D  grows  roughly  as  m  for  large  m ,  with  the  same  precaution  as  in  the  second  derivative 
case.  On  the  other  hand,  the  spectral  radius  of  Dc  grows  as  N2 ,  which  is  again  a  well-known 
property  (see,  for  instance,  [19,  64,  68]).  As  an  example,  Figures  12  and  13  illustrate  for  the 
results  in  Table  8  the  0{c)  and  0(N2)  growth  of  the  spectral  radii  of  D  and  Dc,  respectively.  The 
difference  in  the  asymptotic  behaviors  of  their  spectra  leads  to  different  stability  requirements  in 
the  solution  of  time-dependent  PDEs,  which  we  will  investigate  in  Section  4.2. 
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Figure  9:  Eigenvalues  of  the  first  derivative  matrix  D  incorporating  the  boundary  condition  u(l)  =  0, 
constructed  by  the  scheme  described  in  Section  3.2  using  c  =  207t  and  e  =  10~13. 
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Figure  10:  Spectra  of  D  constructed  with  e  =  10  13  and  (from  (a)  to  (d))  c  =  10,  50, 100, 500. 
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Table  8:  Spectral  radii  p  of  D  constructed  with  £  =  10  13  and  varying  values  of  c.  For  each  c,  the 
spectral  radius  is  compared  with  that  of  the  Chebyshev  matrix  Dq  that  gives  comparable  error  Ec. 


c 

m 

PSWFs 

Ec 

P 

N 

Chebyshe 

Ec 

V 

p 

5 

18 

2.07E-12 

2.63E+01 

23 

7.92E-13 

4.37E+01 

10 

24 

1.27E-12 

4.38E+01 

31 

3.23E-12 

8.06E+01 

20 

34 

4.45E-13 

7.92E+01 

47 

4.37E-13 

1.88E+02 

40 

50 

7.15E-13 

1.43E+02 

73 

5.61E-13 

4.60E+02 

50 

59 

3.91E-13 

1.79E+02 

86 

2.60E-13 

6.41E+02 

80 

81 

2.88E-13 

2.80E+02 

121 

4.74E-13 

1.28E+03 

100 

95 

2.82E-13 

3.45E+02 

144 

4.51E-13 

1.81E+03 

200 

161 

1.93E-12 

6.23E+02 

253 

1.56E-12 

5.63E+03 

400 

292 

1.06E-12 

1.28E+03 

467 

9.41E-13 

1.92E+04 

500 

360 

3.93E-13 

1.64E+03 

575 

4.00E-13 

2.92E+04 

800 

554 

3.43E-13 

2.62E+03 

890 

5.17E-13 

7.00E+04 

1000 

682 

6.23E-13 

3.21E+03 

1090 

1.32E-12 

1.05E+05 

1200 

811 

8.98E-13 

3.81E+03 

1300 

1.05E-12 

1.50E+05 

1600 

1067 

1.14E-12 

5.02E+03 

1710 

1.30E-12 

2.59E+05 

2000 

1324 

1.05E-12 

6.27E+03 

2120 

1.52E-12 

3.98E+05 

N 


Figure  13:  Log-log  plot  of  the  spectral  radius  p  of  the  Chebyshev  differentiation  matrix  Dc  shown  in 
Table  8  against  N. 
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4.2  Numerical  examples 

In  this  subsection,  we  present  the  results  of  several  numerical  experiments  performed  to  assess  the  per¬ 
formance  of  the  scheme  described  in  Section  3.2  in  the  solution  of  time-dependent  partial  differential 
equations,  and  the  associated  eigenvalue  problems.  For  the  solution  of  time-dependent  PDEs,  we  dis¬ 
cretize  the  spatial  derivatives  using  the  differentiation  scheme  described  in  Section  3.2,  and  solve  the 
resulting  system  in  the  time  direction  with  the  explicit  predictor-corrector  schemes  PCI  and  PC3  de¬ 
scribed  in  [26],  where  PCI  is  a  22-step  solver  that  produces  7-digit  precision  at  22  steps  per  wavelength, 
and  PC3  is  a  42-step  solver  that  produces  15-digit  precision  at  36  steps  per  wavelength.  The  predictor- 
corrector  schemes  PCI  and  PC3  are  implemented  with  two  correction  steps  per  time  step,  where  the 
time  steps  are  obtained  from  an  equidistant  discretization  of  the  given  time  interval.  The  starting  values 
for  the  schemes  PCI  and  PC3  are  computed  via  the  spectral  deferred  correction  schemes  DC1  and  DC3 
of  [26],  respectively.  We  refer  the  readers  to  [26]  for  a  detailed  discussion  of  the  construction  of  the 
explicit  schemes  DC1,  DC3,  PCI,  and  PC3,  as  well  as  their  accuracy  and  stability  properties. 

As  the  benchmark,  we  compare  the  differentiation  scheme  described  in  Section  3.2  with  the  Clrebyshev 
collocation  method  and  the  fourth-order  finite  difference  scheme  (FD4).  We  assess  the  performance  of 
both  methods  in  the  solution  of  time-dependent  PDEs  by  combining  them  with  the  time-marching 
schemes  PCI,  PC3,  and  also  the  fourth-order  Runga-Kutta  scheme  (RK4). 

All  schemes  were  implemented  in  FORTRAN  77  and  compiled  with  the  Lahey-Fujitsu  FORTRAN 
95  compiler.  All  experiments  were  run  on  a  PC  with  an  Intel  Core  i7  2.67  GHz  processor  and  12GB 
of  memory,  with  only  one  core  utilized.  No  attempt  was  made  to  parallelize  any  of  the  code.  All 
computations  were  performed  in  double  (16-digit)  precision. 

4.2.1  The  wave  equation 

In  this  subsection,  we  consider  the  following  wave  equation  with  homogeneous  Dirichlet  boundary  con¬ 
dition: 


Utt  —  (4*10) 

u(x,  0)  =  sin(99.57ra;),  Ut(x,  0)  =  0,  (4.11) 

u(0,t)  =  u(2,t)  =  0,  (4.12) 

where  x  £  [0,2]  and  t  >  0.  The  analytical  solution  of  (4.10)-(4.12)  is 

u(x,  t)  =  sin(99.57ra:)  cos(99.57rf).  (4.13) 


By  introducing  the  auxiliary  variable  v  =  Ut,  we  convert  (4.10)  to  a  system  of  two  coupled  equations 
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We  then  solve  the  system  (4.14),  (4.11)-(4.12)  numerically  on  the  time  interval  [0,2.5]  by  combining 
the  differentiation  scheme  described  in  Section  3.2  with  the  time-marching  schemes  PCI  and  PC3.  We 
label  the  combined  schemes  Prolate+PCl  and  Prolate+PC3  respectively.  In  both  of  these  schemes,  the 
differentiation  matrix  was  constructed  with  c  =  315,  which  is  slightly  higher  than  99.57T.  For  compar¬ 
ison,  we  also  solve  the  system  (4.14),  (4.11)-(4.12)  using  a  combination  of  the  Chebyshev  collocation 
method  with  the  PCI,  PC3,  and  RK4  schemes  (which  we  label  Cliebyshev+PCl,  Chebyshev+PC3, 
and  Chebyshev+RK4  respectively),  and  a  combination  of  the  FD4  and  RK4  schemes  (which  we  label 
FD4+RK4). 


Remark  4.3.  In  all  of  the  above  combined  schemes,  we  discretize  the  operator  on  the  right-hand  side  of 
(4.14)  by  the  matrix 


0  I 
D  0 


(4.15) 


where  D  is  the  second  derivative  matrix  incorporating  the  boundary  condition  u(0)  =  u( 2)  =  0,  con¬ 
structed  using  the  corresponding  differentiation  schemes,  and  I  is  an  identity  matrix  with  the  same 
dimension  as  that  of  D.  For  both  the  Chebyshev  collocation  method  and  the  FD4  scheme,  we  incor¬ 
porate  the  boundary  condition  u( 0)  =  u( 2)  =  0  by  first  constructing  the  second  derivative  matrix  as 
described  in  Section  2.3,  and  then  removing  its  first  and  last  rows  and  columns.  Clearly,  in  the  actual 
implementations,  we  only  construct  D ,  and  never  construct  or  apply  (4.15)  explicitly. 


Tables  9-14  summarize  the  results  for  the  schemes  Prolate+PCl,  Prolate+PC3,  Chebyshev+PCl, 
Cliebyshev+PC3,  Chebyshev+RK4,  and  FD4+RK4  respectively.  In  Tables  9  and  10,  c  and  £  denote 
the  bandlimit  and  precision  parameters  used  in  the  differentiation  scheme  described  in  Section  3.2. 
In  all  of  the  tables,  m  denotes  the  number  of  nodes  on  the  interval  [0,  2]  used  in  the  construction  of 
the  differentiation  matrix  D.  In  particular,  it  corresponds  to  D  of  dimension  m  x  m  for  the  scheme 
described  in  Section  3.2,  and  to  D  of  dimension  {m  —  2)  x  (m  —  2)  for  the  Chebyshev  collocation 
method  and  the  FD4  scheme,  p  denotes  the  spectral  radius  of  D;  ns  denotes  the  number  of  nodes  in 
the  equidistant  discretization  of  the  time  interval  [0,2.5],  while  h  denotes  the  corresponding  step-size; 
Ei2  denotes  the  relative  I2  error  of  the  numerical  solution  of  u  at  the  final  time  t  f  =  2.5,  computed  at 
the  spatial  discretization  nodes  on  [0,2],  Finally,  Tc  denotes  the  CPU  times  (in  seconds)  taken  by  the 
scheme  described  in  Section  3.2  to  construct  the  differentiation  matrix,  and  T  denotes  the  CPU  time  (in 
seconds)  taken  by  the  time-marching  scheme. 

In  addition,  we  show  in  Figures  14  and  15  a  comparison  of  the  performance  of  the  above  schemes 
in  the  solution  of  the  system  (4.14),  (4.11)-(4.12).  Figure  14  shows  the  CPU  time  T  versus  the  relative 
error  Ei21  and  Figure  15  shows  the  number  of  time  steps  ns  versus  the  relative  error  E\2.  All  data  points 
in  Figures  14  and  15  are  selected  from  the  results  in  Tables  9-14. 

In  the  following,  we  further  explain  our  results,  and  make  several  observations  and  comments. 


1.  For  each  set  of  results  in  the  tables,  we  first  fix  the  parameters  for  the  differentiation  scheme,  and 
then  increase  the  number  of  time  steps  ns  until  stability  is  attained  and  that  all  observed  error  is 
due  to  the  differentiation  scheme. 
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2.  For  the  scheme  FD4+RK4,  we  do  not  explicitly  construct  the  differentiation  matrix  D  in  (4.15). 
Rather,  we  only  store  the  differentiation  formula  (one  for  differentiation  at  the  first  and  last 
interior  node  on  [0,2],  and  one  for  all  other  interior  nodes)  and  apply  them  one  by  one  during 
actual  computations.  As  a  result,  applying  the  operator  (4.15)  involves  only  0(m)  operations. 
The  differentiation  matrices  constructed  by  Chebysliev  collocation  and  the  scheme  described  in 
Section  3.2  were,  on  the  other  hand,  applied  by  brute-force  multiplications,  and  no  attempt  was 
made  to  speed  up  their  applications. 

3.  For  the  scheme  FD4,  the  differentiation  matrix  D  in  (4.15)  has  real  and  negative  eigenvalues,  and 
its  spectral  radius  grows  as  m2  (see,  for  instance,  [8,  15]).  The  spectral  radii  p  shown  in  Table  14 
are  in  brackets  because  they  are  estimated  by  extrapolation. 

4.  For  the  combined  schemes  that  involve  either  PCI  and  PC3  (see  Tables  9-12),  the  number  of  steps 
ns  is  always  dominated  by  the  stability  requirement.  In  other  words,  once  stability  is  achieved, 
increasing  ns  does  not  increase  accuracy.  From  [26],  the  stability  regions  of  the  schemes  PCI  and 
PC3  are  subsets  of  the  closed  left  half-plane.  We  note  that  the  eigenvalues  of  the  matrix  (4.15)  are 
purely  imaginary,  because  they  are  complex  square  roots  of  the  eigenvalues  of  the  differentiation 
matrix  D,  and  the  latter  eigenvalues  are  all  real  and  negative  (see  Section  4.1.2).  Therefore,  the 
stability  requirement  on  the  step-size  h  is  easily  determined  from  the  span  of  the  stability  regions 
over  the  imaginary  axis  and  the  spectral  radii  p  of  D  shown  in  the  tables.  In  particular,  the  scheme 
PCI  becomes  stable  when  A h  <  0.55,  while  the  scheme  PC3  becomes  stable  when  A h  <  0.33,  where 
A  is  the  maximum  imaginary  part  of  the  eigenvalues  of  the  matrix  (4.15). 

5.  When  compared  to  the  schemes  Chebyshev+PCl,  Chebyshev+PC3,  and 

Chebyshev+RK4,  the  schemes  Prolate+PCl  and  Prolate+PC3  consistently  achieve  speed-up  across 
all  accuracy  levels.  The  speed-up  comes  from  smaller  dimensions  of  differentiation  matrices,  and 
from  smaller  number  of  steps  ns  needed  to  maintain  stability.  For  instance,  from  Tables  10  and 
12,  the  CPU  time  T  it  takes  for  the  scheme  Prolate+PC3  to  achieve  an  accuracy  of  10~12  is  about 
36  times  less  than  it  takes  for  the  scheme  Chebyshev+PC3.  If  the  CPU  time  Tc  taken  in  the 
construction  of  the  differentiation  matrix  is  included,  the  scheme  Prolate+PC3  still  demonstrates 
a  speed-up  of  28  times  compared  to  the  scheme  Chebyshev+PC3. 

6.  The  time  interval  [0,2.5]  contains  about  125  periods  of  the  solution  (4.13),  which  means  it  takes 
about  40  and  120  steps  per  period  for  the  schemes  Prolate+PCl  and  Prolate+PC3  to  achieve  an 
accuracy  of  10-6  and  10  12 ,  respectively.  Again,  these  numbers  of  steps  per  periods  are  higher 
than  those  specified  by  the  schemes  PCI  and  PC3  (see  [26])  because  of  the  stability  restriction 
imposed  by  the  differentiation  scheme. 

7.  From  Tables  9  and  10,  the  CPU  times  Tc  taken  by  the  scheme  of  Section  3.2  to  construct  the 
differentiation  matrices  vary  little  with  the  precision  e,  and  they  are  comparable  in  magnitude  to 
the  marching  times  T.  It  should  be  kept  in  mind,  however,  that  the  times  Tc  are  fixed  given  c  and 
e,  i.e. ,  they  are  independent  on  the  length  of  the  time  interval  tf  on  which  marching  is  performed. 
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The  main  components  of  Tc  come  from  the  construction  of  the  interpolatory  quadrature  using  the 
algorithm  of  [69],  and  from  the  evaluation  of  the  PSWFs.  While  there  are  ways  to  reduce  the  cost 
of  both  components  (in  terms  of  asymptotic  CPU  requirements  and  associated  proportionality 
constants),  we  have  made  no  effort  to  do  so  in  our  implementations. 


Table  9:  Performance  of  the  scheme  Prolate+PCl  in  the  solution  of  (4.14),  (4.11)-(4.12). 


C  £ 

m 

P 

ns 

h 

Eh 

Tc 

T 

315  10_t> 

219 

7.90E+05 

4050 

6.01E-04 

1.04E-05 

1.14 

0.88 

10"6 

222 

9.94E+05 

4550 

5.50E-04 

2.39E-06 

1.14 

0.98 

io-7 

223 

1.23E+06 

5050 

4.95E-04 

9.13E-07 

1.12 

1.12 

10~8 

225 

1.50E+06 

6060 

4.13E-04 

1.60E-06 

1.17 

1.26 

Table  10:  Performance  of  the  scheme  Prolate+PC3  in  the  solution  of  (4.14),  (4.11)-(4.12). 


C  £ 

m 

P 

ns 

h 

Ei2 

Tc 

T 

~3l5  ICT^- 

219 

7.90E+05 

6750 

3.70E-04 

1.04E-05 

1.14 

1.72 

10-6 

222 

9.94E+05 

7550 

3.36E-04 

2.33E-06 

1.14 

1.94 

IO"7 

223 

1.23E+06 

8400 

2.98E-04 

6.26E-07 

1.12 

2.20 

10"8 

225 

1.50E+06 

9300 

2.69E-04 

7.24E-08 

1.17 

2.50 

10"9 

228 

2.14E+06 

11  100 

2.25E-04 

4.28E-10 

1.18 

3.01 

10-10 

230 

2.51E+06 

12  000 

2.08E-04 

3.07E-10 

1.20 

3.24 

10-11 

232 

2.92E+06 

13  000 

1.92E-04 

6.18E-11 

1.22 

3.71 

10-12 

234 

3.37E+06 

14  000 

1.79E-04 

9.53E-12 

1.26 

4.01 

10-13 

238 

3.86E+06 

15  000 

1.67E-04 

1.73E-12 

1.29 

4.31 

Table  11:  Performance  of  the  scheme  Chebyshev+PCl  in  the  solution  of  (4.14),  (4.11)-(4.12). 


m 

P 

ns 

h 

Ei2 

T 

330 

5.55E+08 

108  000 

2.31E-05 

1.28E-03 

52.03 

330 

5.89E+08 

111  000 

2.25E-05 

7.43E-05 

55.42 

340 

6.26E+08 

114  000 

2.19E-05 

2.61E-05 

57.67 

345 

6.63E+08 

117  500 

2.13E-05 

2.48E-05 

62.04 
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Figure  14:  CPU  time  T  versus  relative  error  Ei2  in  the  solution  of  (4.14),  (4.11)-(4.12). 
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Figure  15:  Number  of  time  steps  ns  versus  relative  error  Ei2  in  the  solution  of  (4.14),  (4.11)-(4.12). 
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Table  12:  Performance  of  the  scheme  Chebyshev+PC3  in  the  solution  of  (4.14),  (4.11)-(4.12). 


m 

P 

ns 

h 

Eh 

T 

330 

5.55E+08 

178  500 

2.31E-05 

1.26E-03 

101.65 

335 

5.89E+08 

184  000 

1.36E-05 

8.24E-05 

108.87 

340 

6.26E+08 

189  500 

1.32E-05 

6.85E-06 

111.36 

345 

6.63E+08 

195  500 

1.28E-05 

5.81E-07 

121.38 

350 

7.03E+08 

201  000 

1.24E-05 

3.31E-08 

123.12 

355 

7.44E+08 

207  000 

1.21E-05 

2.41E-09 

135.98 

360 

7.87E+08 

212  500 

1.18E-05 

1.68E-10 

137.81 

365 

8.32E+08 

218  500 

1.14E-05 

1.02E-11 

149.42 

368 

8.59E+08 

222  100 

1.13E-05 

1.27E-12 

156.92 

Table  13:  Performance  of  the  scheme  Chebyshev+RK4  in  the  solution  of  (4.14),  (4.11)-(4.12). 


m 

P 

ns 

h 

Ei2 

T 

330 

5.55E+08 

21  100 

1.12E-04 

1.49E-03 

10.99 

335 

5.89E+08 

22  000 

1.14E-04 

8.79E-05 

11.96 

340 

6.26E+08 

22  400 

1.11E-04 

1.12E-05 

12.34 

345 

6.63E+08 

23  000 

1.09E-05 

8.96E-06 

13.32 

46  000 

5.43E-06 

8.10E-07 

26.94 

350 

7.03E+08 

24  000 

1.04E-04 

7.52E-06 

14.05 

48  000 

5.21E-05 

4.65E-07 

28.18 

60  000 

4.17E-05 

1.92E-07 

35.06 

355 

7.44E+08 

70  000 

3.57E-05 

1.02E-07 

42.80 

90  000 

2.78E-05 

3.76E-08 

55.01 

110  000 

2.27E-05 

1.71E-08 

66.97 

360 

7.87E+08 

80  000 

3.13E-05 

5.98E-08 

49.65 

100  000 

2.50E-05 

2.45E-08 

61.49 

120  000 

2.08E-05 

1.18E-08 

73.70 

150  000 

1.67E-05 

4.85E-09 

92.08 

365 

8.32E+08 

200  000 

1.25E-05 

1.52E-09 

127.69 

250  000 

1.00E-05 

6.22E-10 

159.59 

300  000 

8.33E-06 

2.99E-10 

193.50 

350  000 

7.11E-06 

1.61E-10 

224.30 

368 

8.59E+08 

300  000 

8.33E-05 

3.01E-10 

205.63 

500  000 

5.00E-06 

3.90E-11 

343.58 
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Table  14:  Performance  of  the  scheme  FD4+RK4  in  the  solution  of  (4.14),  (4.11)-(4.12). 


m 

P 

ns 

h 

Eh 

T 

8000 

(1.14E+08) 

10  000 

2.25E-04 

4.20E-04 

5.85 

20  000 

1.25E-04 

1.78E-04 

11.64 

16  000 

(4.54E+08) 

20  000 

1.25E-04 

2.58E-05 

23.68 

40  000 

6.25E-04 

1.11E-05 

47.24 

32  000 

(1.81E+09) 

40  000 

6.25E-04 

1.60E-06 

97.43 

80  000 

3.13E-04 

6.92E-07 

200.82 

64  000 

(7.27E+09) 

80  000 

3.13E-04 

9.74E-08 

435.42 

160  000 

1.56E-04 

4.13E-08 

879.46 
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In  the  rest  of  this  subsection,  we  consider  the  following  wave  equation  with  homogeneous  Dirichlet 
boundary  condition: 


'U'tt  —  ^ xx  i 


(4.16) 


u(x,  0)  =  sin 


,ut(x,0)  =  0, 


(4.17) 


u(  0,  t )  =  u(  2,  t)  =  0, 


(4.18) 


where  x  £  [0,  2],  t  >  0,  and  a  is  a  positive  integer,  of  which  analytical  solution  is 


u(x,  t)  =  sin  (  (  a  +  ^  )  irx  )  cos  (  ( a  +  ^  )  wt  )  . 


(4.19) 


The  value  a  +  1/2  is  the  wavenumber  of  the  solution  u  on  the  spatial  interval  [0, 2]. 

We  present  in  Table  15  the  performance  of  the  scheme  Prolate+PC3  in  the  solution  of  (4.16)-(4.18) 
for  a  range  of  values  of  a.  For  each  a,  we  solve  (4.16)-(4.18)  on  the  time  interval  [0,  tf],  where  tf  is  chosen 
such  that  the  analytical  solution  (4.19)  spans  about  125  periods  on  [0,  £/],  and  compute  the  relative  I2 
error  Ei2  of  the  numerical  solution  at  time  tf.  The  parameters  for  the  scheme  Prolate+PC3  are  chosen 
such  that  about  12  digits  of  accuracy  is  obtained.  Tables  16  and  17  list  for  the  schemes  Chebyshev+PC3 
and  Chebyshev+RK4  the  corresponding  parameters  and  timings  to  obtain  similar  accuracies.  The 
notation  in  Tables  15-17  are  the  same  as  that  in  Tables  9-14. 

In  the  following,  we  make  several  observations  and  remarks  regarding  the  results. 

1.  Again,  for  the  schemes  Prolate+PC3  and  Chebyshev+PC3,  the  number  of  time  steps  ns  is  domi¬ 
nated  by  the  stability  requirement.  For  a  fixed  precision  e,  the  differentiation  matrix  D  in  (4.15) 
constructed  by  the  scheme  described  in  Section  3.2  has  spectral  radius  of  the  order  c2.  As  a  result, 
the  spectral  radius  of  the  matrix  (4.15)  grows  like  c.  Therefore,  for  the  scheme  Prolate+PC3,  the 
time-step  h  decreases  linearly  with  a.  On  the  other  hand,  for  the  scheme  Chebyshev+PC3,  the 
time-step  h  decreases  roughly  as  a2,  due  to  the  0(m4)  growth  of  the  spectral  radius  of  D. 

2.  Compared  to  the  schemes  Chebyshev+PC3  and  Chebyshev+RKT,  the  scheme  Prolate+PC3  demon¬ 
strates  speed-up  across  all  values  of  a.  In  addition,  the  speed-up  increases  with  a  because  of  the 
increasingly  strict  stability  requirement  imposed  by  the  Cliebyshev  differentiation  matrices  as  the 
wavenumber  of  the  problem  increases.  For  a  =  4,  the  CPU  time  T  it  takes  for  the  scheme  Pro- 
late+PC3  to  achieve  an  accuracy  of  10-12  is  about  2.4  times  less  than  it  takes  for  the  scheme 
Chebyshev+PC3;  while  for  a  =  399,  the  speed-up  becomes  about  130  times. 

3.  It  should  be  noted,  however,  that  in  the  computing  environment  in  which  the  experiments  are  run 
(as  described  at  the  beginning  of  Section  4.2),  the  CPU  time  it  takes  for  the  direct  application  of 
an  TV  x  TV  differentiation  matrix  does  not  scale  quadratically  across  all  values  of  N.  In  particular, 
the  proportionality  constant  for  the  direct  application  of  matrices  of  size  N  >  600  is  somewhat 
larger  than  that  associated  with  N  <  600.  This  phenomenon  is  related  to  the  memory  management 
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of  the  particular  computing  environment.  Consequently,  for  a  =  199,  the  scheme  Prolate+PC3, 
with  differentiation  matrix  D  of  size  m  =  445,  demonstrates  a  disproportionate  amount  of  speed¬ 
up  (about  300  times)  when  compared  to  the  scheme  Chebyshev+PC3,  with  D  of  size  N  =  698. 
Therefore,  the  performance  comparison  should  be  seen  more  as  an  indication  than  as  an  absolute 
measure  of  their  relative  performance. 
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Table  16:  Performance  of  the  scheme  Chebyshev+PC3  in  the  solution  of  (4.16)-(4.18) 


a 

tf 

m 

P 

ns 

h 

Ei2 

T 

4 

55.278 

37 

7.99E+04 

47  500 

1.16E-03 

1.49E-12 

0.63 

9 

26.184 

58 

5.01E+05 

56  200 

4.66E-04 

1.67E-12 

1.44 

19 

12.756 

95 

3.70E+06 

74  400 

1.71E-04 

1.71E-12 

4.04 

49 

5.025 

206 

8.37E+07 

131  500 

3.82E-05 

1.35E-12 

26.69 

99 

2.500 

368 

8.59E+08 

222  100 

1.13E-05 

1.27E-12 

156.92 

199 

1.247 

700 

1.13E+10 

402  000 

3.10E-06 

7.42E-13 

4526.4 

399 

0.623 

1350 

1.57E+11 

747  500 

8.33E-07 

1.35E-12 

33  003 

Table  17:  Performance  of  the  scheme  Chebyshev+RK4  in  the  solution  of  (4.16)-(4.18) 


a 

tf 

m 

P 

ns 

h 

Ei2 

T 

4 

55.278 

36 

7.14E+04 

600  000 

9.21E-05 

4.04E-12 

3.03 

9 

26.184 

58 

5.01E+05 

800  000 

3.27E-05 

4.16E-12 

9.78 

19 

12.756 

96 

3.86E+06 

800  000 

1.60E-05 

4.94E-12 

26.83 

49 

5.025 

206 

8.37E+07 

800  000 

6.28E-06 

5.85E-12 

157.68 

99 

2.500 

368 

8.59E+08 

800  000 

3.13E-06 

6.15E-12 

549.76 

199 

1.247 

700 

1.13E+10 

800  000 

1.56E-06 

5.81E-12 

10  953 

399 

0.623 

1350 

1.57E+11 

800  000 

7.79E-07 

6.28E-12 

45  637 

4.2.2  Eigenvalues  of  the  Bessel’s  equation 

In  this  subsection,  we  consider  the  eigenvalues  of  the  Bessel’s  equation: 

1  l 2 

urr  4 — ur - =  Xu,  r  £  [0, 1],  1  =  0,1,2,..., 

with  the  boundary  conditions 


u'(0)  =  0,  w(l)  =  0  if  1  =  0, 
u(0)  =  0,u(l)  =0  if  1^0. 


(4.20) 


(4.21) 

(4.22) 


The  equation  (4.20)  arises  from  the  solution  of  the  wave  equation  in  polar  coordinates.  For  each  l  >  0, 
the  exact  eigenvalues  fe  =  1,2,...  of  (4.20)  are  the  negative  of  the  squares  of  the  roots  of  the  Zt.h 
order  Bessel  function  J/;  in  other  words, 

A  I,k  =  - ilk ,  (4.23) 

where 

Ji(Ji,k)  =  0,  As  =  1,2, -  (4.24) 

The  corresponding  eigenfunctions  are 


Uk{r)  =  Ji{ji,kr),  k  =  1,2, -  (4.25) 

In  the  following,  we  compute  numerically  the  eigenvalues  of  the  Bessel’s  equation  (4.20)  using  the 
scheme  described  in  Section  3.2,  the  Chebyshev  collocation  method,  and  the  fourth-order  finite  difference 
scheme  (FD4),  and  compare  them  to  the  exact  eigenvalues  (4.23).  More  precisely,  we  compute  the 
eigenvalues  of  the  to  x  to  matrix 

M  =  D(2)  +  AD(1)  -  B,  (4.26) 

where  D W  and  D ^  are  the  first  and  second  derivative  matrices  incorporating  the  boundary  conditions 
(4.21)  or  (4.22),  and  A  and  B  are  diagonal  matrices  with  diagonal  entries 

Ai^i  —  1  /  Xi, 

Bi^i  —  l  j Xi,  i  —  1,2,...,  TO 

where  x\, . . . ,  xm  are  the  nodes  used  in  the  construction  of  D W  and  D^. 

Remark  4.4.  For  the  scheme  described  in  Section  3.2,  the  boundary  conditions  (4.21)  for  the  case 
l  =  0  is  incorporated  into  the  orthonormal  set  of  functions  (3.18)  by  subtracting  from  the  PSWFs 
%l>o, ... ,  ipn- 1  the  linear  functions 


(4.27) 

(4.28) 


Hj(x)  =  +  ^'(-i)^  -  !),  j  =  0,  ...,n— 1, 


and  then  applying  the  Gram-Schmidt  algorithm  to  the  matrix  (3.17). 


(4.29) 
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Remark  4.5.  For  the  Chebyshev  collocation  method,  we  incorporate  the  boundary  conditions  (4.21) 
into  D W  and  D ^  of  (4.26)  explicitly  using  the  last  row  and  column  of  the  m  x  m  Chebyshev  matrix 
D  :=  D^'1  given  by  formulas  (2.47)  and  (2.50).  In  particular,  D ^  is  constructed  by  the  formula 

£>(1)  =  Ds  +  uvT,  (4.30) 

where  Ds  is  the  (m  —  2)  x  (m  —  2)  matrix  obtained  by  removing  from  D  its  first  and  last  rows  and 
columns,  and  u  and  v  are  vectors  with  entries: 

fa-i  —  -Ch-f  pm,  (4-31) 

(4.32) 

for  i  =  1, . . . ,  m  —  2.  A  similar  procedure  is  used  to  incorporate  (4.21)  into  the  differentiation  matrices 
of  the  FD4  scheme. 

First,  we  look  at  the  eigenvalues  of  (4.20)  for  the  case  l  =  0.  The  graph  of  the  Bessel  function  Jo 
on  the  interval  [0, 160]  is  shown  in  Figure  16.  Table  18  shows  the  non-zero  eigenvalues  of  the  matrix  M 
in  (4.26)  constructed  using  the  scheme  described  in  Section  3.2  with  c  =  167t  and  s  =  10~13.  There  are 
56  non-zero  eigenvalues,  which  is  equal  to  the  number  of  functions  in  the  orthonormal  set  {</>,;}  (3.18) 
corresponding  to  the  chosen  c  and  e.  The  eigenvalues  Ai, . . . ,  A56  are  real,  distinct,  and  negative;  and 
Ai, . . . ,  A32  are  accurate  to  13  digits  with  respect  to  the  exact  eigenvalues  (4.23).  This  is  expected,  given 
the  choice  c  =  167t  and  the  form  of  the  exact  eigenfunctions: 

Uk{r)  =  Jo(jo.fcr),  (4.33) 

which  contains  k/2  wavelengths  on  the  interval  [0, 1],  for  k  =  1, 2, ... .  For  k  >  32,  the  magnitudes  of  A& 
starts  to  grow  exponentially,  and  they  no  longer  approximate  (4.23). 

For  fixed  e,  the  spectral  radius  p  of  the  matrix  M  grows  as  c2.  As  an  illustration,  Figure  17  plots  the 
spectral  radii  of  M  against  c,  for  £  =  10-7  and  10”13.  All  eigenvalues  of  M  computed  are  either  zero,  or 
real  and  negative.  On  the  other  hand,  the  eigenvalues  of  M  constructed  by  the  Chebyshev  collocation 
method  are  all  real  and  negative.  The  spectral  radii  grow  as  m4,  which  is  illustrated  in  Figure  18. 

Remark  4.6.  For  the  scheme  described  in  Section  3.2,  since  the  prolate  quadrature  nodes  constructed 
with  the  algorithm  of  [69]  do  not  include  the  endpoints  of  the  interval  (see  Remark  3.1),  there  is  no 
division  by  zero  in  the  construction  of  the  matrices  A  and  B  in  (4.26).  On  the  other  hand,  for  the 
Chebyshev  collocation  and  the  FD4  schemes,  we  explicitly  omit  the  endpoints  in  the  construction  of 
D W  and  (see  Remark  4.5),  so  division  by  zero  in  (4.26)  is  avoided. 
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Table  18:  Non-zero  eigenvalues  Ai,...,A56  of  the  matrix  M  constructed  by  the  scheme  described  in 
Section  3.2  using  c  =  167t  and  e  =  10~13.  Ai, . . A32  are  accurate  to  13  digits  with  respect  to  (4.23). 


-.5783185962945060E+01 
-.7488700679069376E+02 
-.2229323036176315E+03 
-.4499335285180342E+03 
-.7558913947839302E+03 
-.1140806031099643E+04 
-.1604677474740231E+04 
-.2147505739697837E+04 
-.2769290832176346E+04 
-.3470032755267547E+04 
-.4249731510652221E+04 
-.5108387099307645E+04 
-.6045999521833543E+04 
-.7062568778614452E+04 
-.8158094869905980E+04 
-.9332577795883790E+04 
-.1058601755667020E+05 
-.1191841415922765E+05 
-.1332976659249901E+05 
-.1482021 199418147E+05 
-.1639970878090875E+05 
-.  1824687473444598E+05 
-.2100242484019345E+05 
-.2578938634166968E+05 
-.3457735371967172E+05 
-.5387757175873770E+05 
-.1040345004919501E+06 
-.4272322530955649E+06 


-.3047126234366095E+02 
-.1390402844264589E+03 
-.3265633529323268E+03 
-.5930428696559549E+03 
-. 9384791 134756928E+03 
-.1362872150854103E+04 
-.1866222004061850E+04 
-.2448528682258052E+04 
-.3109792189768252E+04 
-.3850012528850579E+04 
-.4669189700777149E+04 
-.5567323706308965E+04 
-.6544414545923878E+04 
-.7600462219933994E+04 
-.8735466728550498E+04 
-.9949428071920061E+04 
- .  1 124234625029806E+05 
-.1261422120220677E+05 
-.1406505654372051E+05 
-.1559529716614024E+05 
-.1723462950822729E+05 
-.1928317458550053E+05 
-.2250212957245788E+05 
-.2810174282886851E+05 
-.3863584374915628E+05 
-.6174489536835697E+05 
-.1276690375260464E+06 
-.4908888139141377E+06 
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Figure  17:  Log-log  plot  of  the  spectral  radius  p  of  the  matrix  M  constructed  by  the  scheme  in  Section  3.2 
against  c,  with  e  =  10~7, 10”13. 

Next,  we  look  at  the  convergence  performance  of  the  scheme  described  in  Section  3.2,  the  Clrebyshev 
collocation  scheme,  and  the  FD4  scheme  in  approximating  the  exact  eigenvalues  (4.23)  of  the  Bessel’s 
equation.  The  column  labeled  ‘PSWFs’  in  Table  19  shows  the  relative  errors  Erei  when  the  eigenvalue 
Ao,4o  is  computed  using  the  scheme  described  in  Section  3.2  with  e  =  10-13  and  varying  values  of 
c.  For  comparison,  Table  19  also  shows  the  errors  Erei  when  Aoqo  is  computed  using  the  Clrebyshev 
collocation  and  the  FD4  schemes.  Also,  Figure  19  shows  the  plot  of  the  errors  Erei  in  Table  19  against 
the  number  of  nodes  to.  In  addition,  Figures  20-22  show  the  plots  of  the  errors  Erei  against  to  when  the 
eigenvalues  Ayioi  Aio.ioo,  and  A50,200  are  computed  using  the  above  three  schemes.  From  Figures  19- 
22,  we  see  that  both  the  scheme  described  in  Section  3.2  and  the  Clrebyshev  collocation  scheme  display 
exponential  (spectral)  convergence  in  accuracy.  However,  considering  the  frequency  of  the  corresponding 
eigenfunctions,  the  scheme  described  in  Section  3.2  requires  fewer  points  per  wavelength  to  achieve  double 
precision  accuracy.  For  instance,  for  the  computation  of  Aso^oco  it  takes  about  2.6  points  per  wavelength 
for  the  scheme  described  in  Section  3.2  to  achieve  15  digits  of  accuracy,  while  the  Chebyshev  collocation 
scheme  takes  about  4  points  per  wavelength  to  achieve  the  same  order  of  accuracy.  The  FD4  scheme, 
on  the  other  hand,  only  displays  1  to  2  digits  of  accuracies  with  comparable  numbers  of  points. 

Remark  4.7.  We  would  like  to  reiterate  that  for  the  scheme  described  in  Section  3.2,  the  number  of  nodes 
to.  is  determined  by  the  choice  of  c  and  e.  We  plot  in  Figures  19-22  the  relative  errors  Erei  obtained  by 
the  scheme  against  m  for  ease  of  comparison  to  the  Chebyshev  collocation  and  the  FD4  schemes. 
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Figure  18:  Log-log  plot  of  the  spectral  radius  p  of  the  matrix  M  constructed  by  the  Chebysliev  collocation 
method  against  the  number  of  nodes  m. 


Table  19:  Relative  errors  Erei  of  computing  Ao,4o  using  the  scheme  described  in  Section  3.2,  the  Cheby- 
shev  collocation  scheme,  and  the  FD4  scheme. 


c 

PSWFs 

TO  Erei 

Chebyshev 
to  Erei 

m 

FD4 

Erei 

48 

56 

1.01E-03 

68 

1.34E-03 

56 

2.02E-01 

50 

59 

3.60E-05 

72 

3.39E-05 

63 

1.51E-01 

52 

60 

8.38E-07 

76 

4.69E-07 

70 

1.17E-01 

54 

62 

1.94E-07 

80 

9.92E-09 

77 

7.66E-02 

56 

63 

9.92E-08 

84 

3.46E-10 

84 

4.64E-02 

58 

65 

2.57E-10 

88 

1.05E-11 

91 

3.29E-02 

60 

66 

1.52E-11 

92 

2.54E-13 

98 

2.42E-02 

62 

68 

5.72E-15 

96 

3.15E-15 

105 

1.82E-02 

57 


PSWFs 
Chebyshev 
‘  FD4 


Erel 


Erel 


X 

X 


Figure  19:  Relative  errors  Erej  against  m  in  the  computation  of  Aq.4o- 
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Figure  20:  Relative  errors  Erei  against  m  in  the  computation  of  Ai.io 


4.2.3  A  variable-coefficient  hyperbolic  equation 

In  this  subsection,  we  consider  the  following  variable-coefficient  hyperbolic  equation: 


ut  =  a(x)ux, 

u(x,  0)  =  sin(307rx  +  1.5  sin(8.97rx)), 
=  g(t), 

where  x  £  [—  1, 1],  t  >  0,  and 


307r  +  13.357T  cos(8.97rx)  ’ 
g(t)  =  sin(307r  +  1.5  sin(8.97r)  +  t). 

The  analytical  solution  of  (4.34)-(4.38)  is  given  by 

u(x,  t)  =  sin(307ra;  +  1.5sin(8.97rx)  +  t). 


(4.34) 

(4.35) 

(4.36) 


(4.37) 

(4.38) 


(4.39) 


Using  the  transformation 


v(x,t)  =  u(x,t)  -  g(t), 


we  reduce  (4.34)-(4.36)  into  a  hyperbolic  equation  with  zero  boundary  condition  at  x  =  1: 


(4.40) 


vt  =  a(x)vx  -  (4.41) 

v(x,  0)  =  sin(307ra:  +  1.5  sin(8.97ra;))  —  g( 0),  (4.42) 

v(l,  t)  =  0.  (4.43) 


In  the  following,  we  solve  (4.41)-(4.43)  numerically  on  the  time  interval  [0,  1000]  using  the  schemes 
Prolate+PCl,  Prolate+PC3,  Clrebyshev+PCl,  Chebyshev+PC3,  and  Chebyshev+RK4,  and  compare 
their  timings  and  accuracies.  For  the  schemes  Prolate+PCl  and  Prolate+PC3,  we  construct  the  matrix 
M  discretizing  the  operator  a  ■  vx  by  first  constructing  an  m  x  rn  first  derivative  matrix  D  incorporating 
the  boundary  condition  u(  1)  =  0  using  the  scheme  described  in  Section  3.2,  and  then  multiplying  D  on 
the  left  by  the  m  x  m  diagonal  matrix  A  with  diagonal  entries 


Aiti  =  a(xi),  i  = 


(4.44) 


where  x\, . . . ,  xm  are  the  nodes  used  in  the  construction  of  D.  The  matrices  M  are  similarly  constructed 
for  the  schemes  Chebyshev+PCl,  Chebyshev+PC3,  and  Chebyshev+RK4,  in  which  we  construct  the 
first  derivative  matrix  Dq  using  the  Chebyshev  collocation  method  described  in  Section  2.3.2,  and  then 
strip  Dc  of  its  first  row  and  column. 

Tables  20-24  summarize  the  results  for  the  schemes  Prolate+PCl,  Prolate+PC3,  Chebyshev+PCl, 
Clrebyshev+PC3,  and  Chebyshev+RK4  respectively.  In  the  tables,  m  denotes  the  number  of  nodes 


60 


on  [—1, 1]  used  in  the  construction  of  the  matrix  M,  p  denotes  the  spectral  radius  of  the  discretized 
operator  a  ■  vx,  and  Ei2  denotes  the  relative  I2  error  of  the  numerical  solution  of  v  at  the  final  time 
tf  =  1000.  The  rest  of  the  notation  is  the  same  as  that  in  Tables  9-13  in  Section  4.2.1.  Figure  23 
shows,  for  the  above  schemes,  the  CPU  time  T  versus  the  relative  error  E/,2 ,  while  Figure  24  shows  the 
number  of  time  steps  ns  versus  Ei2.  All  data  points  in  Figures  23  and  24  are  selected  from  the  results 
in  Tables  20-24.  In  addition,  Figures  25  and  26  show,  for  selected  sets  of  parameters,  the  spectra  of  the 
matrix  M  constructed  using  the  scheme  described  in  Section  3.2  and  the  Chebyshev  collocation  method, 
respectively. 

In  the  following,  we  make  several  observations  and  comments  based  on  the  presented  results,  and  on 
the  results  of  our  more  extensive  experiments. 

1.  In  all  of  the  above  schemes,  the  eigenvalues  of  the  matrix  M  lie  on  the  left  half-plane,  so  stability 

is  guaranteed  provided  that  sufficiently  small  time-steps  are  chosen  for  the  time-marching  schemes 
PCI,  PC3,  and  RK4.  In  particular,  for  the  schemes  Prolate+PCl,  Prolate+PC3,  Chebyshev+PCl, 
and  Chebyshev+PC3,  the  number  of  steps  ns  is  always  dominated  by  the  stability  requirement. 
On  the  other  hand,  for  the  scheme  Chebyshev+RK4,  ns  is  dominated  by  the  stability  requirement 
up  to  a  desired  accuracy  of  about  10  ,  after  that,  ns  is  dominated  by  the  accuracy  requirement. 

2.  Compared  to  the  schemes  Chebyshev+PCl,  Chebyshev+PC3,  and  Chebyshev  +RK4,  the  schemes 
Prolate+PCl  and  Prolate+PC3  are  superior  across  all  accuracies,  both  in  terms  of  CPU  time  T 
and  the  number  of  time  steps  ns.  This  is  because,  as  indicated  in  Tables  20-24,  the  matrices  M 
constructed  by  the  scheme  described  in  Section  3.2  have  smaller  spectral  radii  p  compared  to  those 
constructed  using  the  Chebyshev  collocation  method,  given  the  same  accuracy  requirement.  In 
particular,  our  more  extensive  experiments  show  that  the  matrices  M  constructed  by  the  former 
scheme  have  p  that  grow  as  c  for  fixed  e,  while  those  constructed  by  the  latter  scheme  have  p  that 
grow  as  m. 

3.  From  Table  21,  we  see  that  the  bandlimit  parameter  c  required  to  solve  for  v  in  (4.40)  to  precision 
10~13  is  about  480.  This  corresponds  to  about  153  wavelengths  on  the  interval  [—1, 1].  Thus,  the 
frequency  of  v  in  the  spatial  dimension  is  much  higher  than  its  frequency  in  the  time  dimension, 
the  latter  of  which  equals  1/27 r.  As  a  result,  the  penalty  imposed  by  the  stability  requirement 
in  this  example  is  higher  than  that  in  the  example  of  Section  4.2.1.  In  particular,  the  solution  v 
spans  about  160  periods  on  the  time  interval  [0,1000],  which  means  that  it  takes  about  130  and 
590  steps  per  period  for  the  schemes  Prolate+PCl  and  Prolate+PC3  to  achieve  an  accuracy  of 
5x  10~'  and  10— 13 ,  respectively.  These  results,  although  less  than  optimal,  still  compare  favorably 
to  those  obtained  by  the  schemes  Chebyshev+PCl,  Chebysliev+PC3,  and  Chebyshev+RK4. 


61 


Prolate+PCl  — Ef 
Prolate+PC3  — Ar 
Chebyshev+PCl  — ©- 
Chebysliev+PC3  — $- 
Chebyshev+RK4  — *- 


104  | 
103  E 

102  f 

101  f 

10°  - 
10' 


10 


-4 


10 


-6 


10"8 

Eu 


10 


-10 


10- 


Figure  23:  CPU  time  T  versus  relative  error  Ei2  in  the  solution  of  (4.41)-(4.43). 
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Figure  24:  Number  of  time  steps  ns  versus  relative  error  Ei2  in  the  solution  of  (4.41)-(4.43). 
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Figure  25:  Spectra  of  M  constructed  using  the  scheme  described  in  Section  3.2,  with  e  =  10  13  and 
(from  (a)  to  (d))  c  =  200,400,600,800. 
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Figure  26:  Spectra  of  M  constructed  using  the  Chebyshev  collocation  scheme,  with  (from  (a)  to  (d)) 
to  =  200,400,600,800. 
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Table  20:  Performance  of  the  scheme  Prolate+PCl  in  the  solution  of  (4.41)-(4.43). 


c 

£ 

m 

P 

ns 

h 

Ei2 

Tc 

T 

300 

10-8 

217 

1.15E+01 

21  000 

4.76E-02 

1.20E-05 

0.79 

4.12 

310 

10"8 

226 

1.16E+01 

21  200 

4.72E-02 

4.61E-07 

0.86 

4.24 

320 

10"8 

230 

1.18E+01 

21  400 

4.67E-02 

3.69E-07 

0.92 

4.52 

330 

10"8 

237 

1.19E+01 

21  600 

4.63E-02 

5.44E-07 

1.01 

4.86 

Table  21:  Performance  of  the  scheme  Prolate+PC3  in  the  solution  of  (4.41)-(4.43) 

c 

£ 

m 

P 

ns 

h 

Ei2 

T 

T 

330 

10"y 

242 

1.41E+01 

43  000 

2.33E-02 

1.61E-06 

1.02 

10.56 

350 

10-9 

252 

1.44E+01 

44  000 

2.27E-02 

7.00E-08 

1.12 

11.53 

370 

10-9 

268 

1.58E+01 

48  000 

2.08E-02 

3.36E-09 

1.32 

14.24 

390 

10-11 

283 

2.00E+01 

61  000 

1.64E-02 

6.28E-10 

1.51 

20.63 

410 

10-11 

299 

2.04E+01 

62  000 

1.61E-02 

3.79E-10 

1.70 

23.64 

420 

10~n 

305 

2.05E+01 

62  500 

1.60E-02 

1.11E-11 

1.76 

24.57 

430 

10-13 

314 

2.51E+01 

76  000 

1.32E-02 

5.15E-12 

1.90 

30.50 

450 

10"13 

330 

2.55E+01 

77  500 

1.29E-02 

2.31E-12 

2.09 

34.55 

470 

10-13 

343 

2.74E+01 

83  500 

1.20E-02 

4.99E-13 

2.31 

41.26 

480 

10-13 

349 

2.76E+01 

84  000 

1.19E-02 

1.30E-13 

2.42 

42.56 

Table  22:  Performance  of  the  scheme  Chebyshev+PCl  in  the  solution  of  (4.41)-(4.43). 


m 

P 

ns 

h 

Eh 

T 

280 

1.27E+02 

226  000 

4.42E-03 

1.06E-05 

72.21 

300 

1.45E+02 

256  000 

3.91E-03 

3.43E-06 

93.32 

320 

1.66E+02 

291  500 

3.43E-03 

3.81E-06 

120.92 

340 

1.87E+02 

330  000 

3.03E-03 

4.16E-06 

154.23 

Table  23:  Performance  of  the  scheme  Chebyshev+PC3  in  the  solution  of  (4.41)-(4.43). 


m 

P 

ns 

h 

Eh 

T 

300 

1.45E+02 

430  000 

2.33E-03 

9.69E-07 

173.22 

320 

1.66E+02 

490  000 

2.04E-03 

3.77E-07 

226.92 

340 

1.87E+02 

560  000 

1.79E-03 

7.28E-07 

289.54 

360 

2.10E+02 

620  000 

1.61E-03 

6.32E-09 

358.53 

380 

2.34E+02 

690  000 

1.45E-03 

9.84E-10 

435.70 

400 

2.59E+02 

770  000 

1.30E-03 

2.73E-10 

542.54 

420 

2.86E+02 

855  000 

1.17E-03 

2.77E-11 

776.64 

440 

3.14E+02 

940  000 

1.06E-03 

2.09E-12 

911.67 

460 

3.43E+02 

1  040  000 

9.62E-04 

5.11E-13 

1014.6 

480 

3.73E+02 

1  130  000 

8.85E-04 

1.52E-13 

1328.2 
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Table  24:  Performance  of  the  scheme  Chebyshev+RK4  in  the  solution  of  (4.41)-(4.43). 


m 

P 

ns 

h 

El  2 

T 

280 

1.27E+02 

45  000 

2.22E-02 

1.00E-05 

17.25 

300 

1.45E+02 

52  000 

1.92E-02 

9.78E-07 

23.81 

320 

1.67E+02 

58  000 

1.72E-02 

3.99E-07 

28.92 

340 

1.87E+02 

66  000 

1.52E-02 

8.56E-08 

37.97 

360 

2.10E+02 

74  000 

1.35E-02 

2.11E-08 

48.39 

380 

2.34E+02 

82  000 

1.22E-02 

1.36E-08 

59.50 

102  000 

9.80E-03 

5.89E-09 

75.02 

400 

2.59E+02 

105  000 

9.52E-03 

4.96E-09 

84.34 

200  000 

5.00E-03 

5.21E-10 

159.42 

420 

2.86E+02 

200  000 

5.00E-03 

3.74E-10 

214.88 

300  000 

3.33E-03 

8.19E-11 

321.41 

400  000 

2.50E-03 

3.87E-11 

429.48 

440 

3.14E+02 

500  000 

2.00E-03 

1.01E-11 

574.66 

750  000 

1.33E-03 

3.11E-12 

867.72 

460 

3.43E+02 

750  000 

1.33E-03 

2.11E-12 

788.76 

1  000  000 

1.00E-03 

9.33E-13 

1054.3 

480 

3.74E+02 

750  000 

1.33E-03 

1.90E-12 

1055.8 

1  500  000 

6.67E-04 

1.62E-13 

2086.1 
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5  Conclusions  and  future  work 


We  have  presented  a  new  class  of  numerical  differentiation  schemes  constructed  via  the  PSWFs.  As 
opposed  to  existing  collocation  methods,  the  schemes  are  based  on  the  construction  of  an  approximate 
interpolation  u  of  a  function  /  via  a  least-squares  type  procedure,  in  which  we  do  not  require  u  to  be 
exactly  equal  to  /  at  any  of  the  interpolation  nodes.  For  problems  that  involve  bandlimited  functions, 
the  schemes  require  fewer  points  per  wavelength  to  attain  the  same  accuracy  when  compared  to  the 
Clrebyshev  collocation  method.  In  addition,  when  solving  time-dependent  PDEs  with  non-periodic 
boundary  conditions,  the  resulting  first  and  second  derivatives  have  spectral  radii  that  grow  as  c  and 
c2  respectively,  for  fixed  precision  e.  Our  numerical  experiments  indicate  that,  when  combined  with  a 
numerical  ODE  solver  to  solve  time-dependent  PDEs,  the  schemes  outperform  the  Chebyshev  collocation 
and  the  finite  difference  methods,  in  particular  when  high  accuracy  is  required  or  the  solutions  contain 
large  numbers  of  wavelengths. 

In  the  following,  we  discuss  several  possible  extensions  to  the  schemes: 

1.  It  is  possible  to  accelerate  the  application  of  the  differentiation  matrices  constructed  by  the  schemes 
to  vectors.  In  particular,  the  differentiation  matrices  D  of  the  schemes  presented  here  satisfy  the 
symmetric  property 

D ,  j  —  Dm_2)Tn_j,  1  A  j,  j  ft  m, 
in  the  case  of  first  derivatives,  and 

D  j  j  —  Dm—irn—ji  1  f:  b  J  — 

in  the  case  of  second  derivatives,  where  in  is  the  dimension  of  D.  Thus,  a  fast  algorithm  for 
applying  I?  to  a  vector  v  that  is  similar  to  the  scheme  of  [60]  can  be  constructed. 

2.  On  the  other  hand,  the  schemes  presented  here  can  be  modified  to  make  the  resulting  differentia¬ 
tion  matrices  amenable  to  fast  applications.  For  pseudospectral  methods  based  on  Chebyshev  or 
Legendre  polynomials,  the  fast  Fourier  transform  (FFT)  or  the  fast  multipole  method  (FMM)  (see, 
for  instance,  [12,  31,  33])  can  be  employed  to  reduce  the  cost  of  applying  an  N  x  N  differentiation 
matrix  to  a  vector  to  0(N  log  N)  operations  (see  [8,  20,  63]).  For  pseudospectral  methods  based  on 
the  PSWFs,  the  authors  of  [39]  constructed  an  algorithm  that  utilizes  the  fast  multipole  method 
(FMM)  to  reduce  the  cost  to  O(N)  operations.  Similar  modifications  to  our  schemes  are  currently 
under  investigation  by  the  authors,  and  any  results  will  be  reported  in  a  later  date. 

3.  For  the  schemes  presented  here,  the  interpolation  nodes  on  which  the  differentiation  matrix  D  is 
constructed  all  lie  in  the  interior  of  the  interval  [—1,1]  (see  Remark  3.1),  and  boundary  condi¬ 
tions  are  incorporated  implicitly  in  the  orthonormal  set  of  functions  (f>i, . . . ,  (j>k  (see  (3.14)).  This 
may  seem  to  pose  difficulty  in  some  two-dimensional  problems,  such  as  the  two-dimensional  wave 
equation 

Utt  —  T  n yy,  1  b  Xf  y  b  1,  (h.l) 
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with  boundary  conditions 


u(x,  y) 


sin(7ra;)  for  y  =  ±1, 
0  for  x  =  ±1, 


(5.2) 


since  then  the  boundary  conditions  cannot  be  incorporated  independently  into  the  differentiation 
matrices  discretizing  uxx  and  uyy. 

One  possible  way  to  tackle  it  is  to  modify  the  differentiation  matrix  D  in  (3.22): 


D  =  UP*W 


(5.3) 


into  the  form 

D  =  YDX,  (5.4) 

where  X  is  a  matrix  that  interpolates  from  a  set  of  points  yi , . . . ,  yi  to  the  quadrature  nodes 
x\, . . . ,  xm  used  in  the  construction  of  D ,  with  containing  the  end  points  of  the  interval 

[—1,1];  and  Y  is  a  matrix  that  interpolates  from  x±, . . . ,  xm  to  y\, . . . ,  yi .  Both  X  and  Y  can  be 
constructed  via  least-squares.  By  using  the  matrix  D  instead  of  D  in  a  time-marching  scheme,  the 
boundary  conditions  can  be  enforced  directly  at  each  time-step. 

In  addition,  using  D  as  the  differentiation  matrix  allows  the  incorporation  of  boundary  conditions 
via  the  method  described  in  Remark  4.5,  which  can  be  more  convenient  in  certain  problems. 
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